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ABSTRACT 


The control of an airbreathing turbojet engine by an on-board 
digital computer is studied. The approach taken is to model the turbo- 
jet engine as a linear, multivariable system whose parameters vary with 
engine operating environment. From this model adaptive closed-loop or 
feedback control laws are designed and applied to the acceleration of 
the turbojet engine. 

A linear state variable model of turbojet engine dynamics is identi- 
fied by a technique that determines first the model structure then the 
model parameters. Models are identified at several operating conditions 
to completely describe the entire engine operating range. Only inputs 
and noise corrupted outputs realizable at an actual engine are con- 
sidered . 

Adaptive feedback controls are designed using sampled-data control 
theory. The necessary optimality conditions for the optimal sampled- 
data output regulator are derived. These necessary conditions and a 
variable sampling rate to reduce computer processing time (adaptive 
sampling) are combined to form an adaptive digital control scheme. 

This scheme generates constant proportional feedback laws that are 
functionally dependent on the control system sampling rate and that 
require process outputs rather than states for control purposes. This 
adaptive digital control scheme is applied to the previously identified 
engine models to control a turbojet engine. Several engine accelera- 
tion transients are simulated to study the effectiveness of the result- 
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ant adaptive control and the relative improvements in computer proc- 
essing time. Additionally, the incorporation of certain physical en- 
gine constraints into the overall control problem is considered. 
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CHAPTER I 


INTRODUCTION AND BACKGROUND 

The purpose of this dissertation Is to develop techniques for the 
application of modern control theory to turbojet engine control system 
design. The approach taken Is to consider the turbojet engine as a 
linear, multivariable, dynamically varying system and design adaptive 
feedback controls that meet engine operation and performance require- 
ments. Specific techniques used in the design process Include stochas- 
tic system Identification, discrete output regulator theory, and adapt- 
ive sampling. Before proceeding with the development and application 
of the adaptive control design, the history and significance of turbo- 
jet engine control are discussed. 

1.1 Engine Control Problem 

Initially turbojet engine configurations were simple combinations 
of a compressor, combustor, turbine, and exhaust nozzle. As Sobey and 
Suggs (1963) indicate, the first control systems for these engines were 
hydromechanical and used the principle of Che flyball governor exten- 
sively for fuel-rotor speed control. As performance demands on turbo- 
jet engines increased, so did engine and control system complexity. 

The addition of a second compressor, driven by its own turbine, enabled 
greater flexibility of compressor performance at high discharge to in- 
take pressure ratios. However, this "twin spool" arrangement put 
greater demands on the control system since the rotational speed of two 
mechanically Independent turboshafts were now to be controlled. The 
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addition of thrust augmentation schemes such as afterburning and vari- 
able exhaust nozzle areas added auxiliary control tasks to the basic 
fuel-speed control. With increased performance demands physical engine 
constraints such as maximum allowable turbine temperature and stable 
compressor operation became important engine control system design con- 
siderations. Still another major factor in the evolution of engine con- 
trol systems was the application of engines to advanced aircraft propul- 
sion. An example would be the application of the variable cycle engine 
to supersonic or short takeoff and landing (STOL) aircraft (Beattie, 
1974). Since a variable cycle engine incorporates variable compressor 
and turbine geometries on both spools and two variable area exhaust noz- 
zles, it requires more control functions than current engines and there- 
fore will require a correspondingly complex control. 

Thus from a control viewpoint a modern engine can be considered as 
a nonlinear multivariable (multi -input multi-output) system with sev- 
eral different control tasks to be accomplished simultaneously. In 
spite of this complexity, however, each control system can be consid- 
ered as consisting of a basic fuel-speed control and a variety of other 
auxiliary control functions. The basic requirements of such a fuel- 
speed control are (1) to accelerate the engine without violation of 
physical constraints and (2) to control steady-state fuel flow. Al- 
though many schemes have been developed to accomplish these basic ob- 
jectives, they all contain the basic structure shown in Figure 1.1. 

The computational device is given as input information the com- 
manded throttle setting, the environmental conditions of the engine 
(e.g,, altitude and flight speed), and some physical engine variables 
(rotor speeds, temperatures, pressures, etc.). From this information 
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Engine outputs 



Figure 1-1. - Basic fuel-speed control system. 
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the computational device geneirates the required fuel metering valve po- 
sition and consequently the appropriate engine fuel flow. The computa- 
tional device is generally mechanized as either a hydromechanical or an 
electrical system, or a combination of both. 

Hydromechanical controls are the oldest and most popular of the 
available mechanizations. As Leeson (1974) points out hydromechanical 
controls are essentially devices that maintain a schedule of desired 
engine temperature or acceleration. The scheduling and resultant multi 
plication can be accomplished in a variety of ways using cams, linkages 
nozzles, springs, and valves. On the other hand steady-state fuel con- 
trol is normally accomplished by simple flyball type governors that me- 
chanically sense rotor speed and directly move the metering valve. 

Electronic fuel controls are thought of as a modern innovation. 
However, one of the first serious applications was developed in the 
late 1940 *s (Leeson, 1974). Electronic controls can be divided into 
two categories, analog and digital. The vast majority of electronic 
engine control systems are analog. Typically an analog control would 
perform the same control functions as its hydromechanical counterpart 
but with electromagnetic pickups for sensing rotor speed, thermocouples 
for sensing temperature , differential transformers for position indica- 
tion, a two-stage servo valve to perform the necessary work for meter- 
ing valve position, and various electronic amplifiers. An example of 
an analog engine control is given in detail by Prue (1974) and Loft 
(1969). 

Althoug^i most controls are either hydromechanical or analog, most 
of the current research interest and emphasis is in digital electronic 
engine control because of its future promise. Digital engine control 
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requires the use of a digital computer with either supervisory control 
over a basic hydromechanical fuel control system (Griffiths and Powell, 
1974) or full authority control over all aspects of engine control func- 
tions (Bentz, 1974). Actual applications of digitally controlled jet 
engines have been reported by Cwynar and Batterton (1975), Batterton 
et al. (1974), Waters (1974), Arpasi et al. (1972), Frazzini (1970), 
Eccles and Shutler (1970), Bayati and Frazzini (1968); and Lewis and 
Munns (1968). The next section discusses how engine control systems, 
both hydromechanical and electronic, have been designed in the past and 
possible future design techniques . 

1.2 Engine Control Design Methods 

For a complex engine configuration much of the control development 
is based on good steady-state turbomachinery operation and acceleration 
response with respect to throttle changes while maintaining various 
physical engine constraints. Traditionally the control requirements 
were met by scheduling engine inputs as functions of flight conditions, 
pilot throttle demand, and one engine output, rotor speed. Recently, 
however, control systems have been designed that use additional meas- 
ured engine output variables to yield better steady-state and transient 
definitions of the engine operating constraints. The vast majority of 
these systems have been designed using classical frequency response 
techniques . 

Classical frequency response techniques are restricted to single- 
input single-output time-invariant systems. Consequently for a multi- 
variable engine the control function for each input is designed inde- 
pendently. When these independently designed control functions are 
combined into a complete engine control, input interaction may slgnifi- 
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cantly degrade engine performance. This problem can be overcome by ad- 
justing appropriate control function bandwidths or by designing decou- 
pling paths between interacting engine inputs . However , bandwidth ad- 
justments degrade system response^ and decoupling, a nonsystematic pro- 
cedure, may require many attempts to find an acceptable solution. In 
addition to these traditional control problems increased engine com- 
plexity and performance requirements (e.g., better integration of all 
engine control functions, optimization of fuel consumption, etc.) have 
placed demands on the development of control systems that cannot be met 
by traditional design methods. Modem Control Theory (MCT) offers pos- 
sible solutions to these design problems . 

MCT is a general title that includes several different control 
concepts. Some of these concepts are the state space representation of 
systems, optimal control theory, estimation and Identification theory, 
Pontryagin*s maximum principle, and several vector frequency response 
tehcniques. In general MCT design techniques are computer oriented and 
thus can systematically handle more complex multivariable problems . 

MCT has therefore become an increasingly important tool to many indus- 
trial and research concerns in the design and analysis of jet engine 
controls . 

In particular, some preliminary engine control design using vector 
frequency response techniques has been done by MacFarlane, et al. 

(1971) and McMorran (1970). Chen (1972), Ahlbeck (1966), and Mueller 
(1971) have applied frequency response techniques to find the transfer- 
function matrix of known dimension of a gas turbine system. Also, 
Michael and Farrar (1973) have applied continuous -time linear optimal 
state regulator theory (Kwakernaak and Sivan, 1972) to engine control 
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design. This work asstjmed the availability of each state to implement 
the control law. Michael and Farrar (1973) have used a least squares 
curve fitting technique to fit an assumed model to engine simulation 
data. Recently, Michael and Farrar (1975) combined their least squares 
identification with a dynamic nonlinear filter to identify gas turbine 
dynamics from stochastic input-output data. In each of the identifica- 
tion papers a priori assumptions were made about system order and 
structure. Sevich and Beattie (1975) have used nonlinear programming 
to develop optimal engine variable trajectories. 

Much of the potential of MCT has not been fully realized. The 
principal objective of this dissertation is to develop design tech- 
niques that further exploit the capabilities of MCT when applied to 
jet engine control. In the next section research areas of signifi- 
cance are identified, problem objectives defined, and the proposed so- 
lutions outlined. 

1.3 Problem Statement 

One of the basic assumptions of this research is the presence of 
an on-board digital computer for full authority engine control. Since 
there exists a finite limit to the time available for control update 
purposes, the efficiency with idiich the computer functions is of the 
utmost importance. Efficient computer utilization would allow (1) the 
time-sharing of several control tasks by a single computer or (2) the 
use of small, less expensive, specialized computers. 

At first, complex continuous -time control systems were approxi- 
mated on the digital computer (as in Michael and Farrar, 1973), Such 
approximation techniques often are computationally inefficient and re- 
quire large, fast, expensive machines to achieve a satisfactory ap- 
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proxlmation. Alternately, since a computer accepts data in discrete 
form, a control designed by sampled-data theory coul^ be used. 

When a digital computer is introduced into the control loop, the 
resultant system can be handled by sampled-data control theory (Kuo, 
1970). Sampled-data theory allows the control designer to implement a 
discrete equivalent of the continuous solution rather than an approxi- 
mation of the continuous solution with a digital computer. A control 
system designed by sampled-data theory would allow (1) the use of 
smaller, less expensive computers and (2) the utilization of computer 
time-sharing capabilities (Levis, et al., 1971) and would therefore 
provide for efficient use of the computational facilities. To further 
increase the efficiency of the control computer, adaptive sampling 
(Dorf, et al., 1962) can be Introduced into the control algorithm. 

Adaptive sampling varies the frequency with which the computer 
samples the continuous signal for digital processing. The frequency 
is varied as a function of some continuous system parameter. The 
overall effect of such a scheme is to Increase the control activity of 
the computer during high Information periods (engine transients, for 
example) and reduce the activity during periods of low information. 
Thus the first objective of this dissertation is to develop a sampled- 
data (discrete) engine control algorithm that incorporates adaptive 
sampling for efficient on-board computer utilization. 

In addition to efficient computer operation there is a need for a 
systematic design procedure that yields a practical and Implementable 
control law and eliminates the problem of input Interaction in a com- 
plex engine. Michael and Farrar (1973) have shown that an adaptive 
control designed by continuous -time state regulator theory can fulfill 
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this requirement. In general regulator theory can be used not only to 
systematically design controls that take advantage of input interac- 
tion, but also to easily evaluate control effectiveness. Also, the 
linear feedback law of optimal regulator theory is both practical and 
implementable . 

One drawback of the control formulation of Michael and Farrar is 
the assumption of full state availability for control purposes. The 
resultant control is Inflexible in that state variables must be physi- 
cally present in either a sensed or estimated form. One alternative 
to this is the linear output regulator formulation of Levine and Athans 
(1970), This output regulator formulation retains the benefits of 
state regulator theory but no longer requires f\xll state availability. 
Different combinations and numbers of output variables can be used as 
feedback variables and the resultant control laws readily designed. 
Thus, the second objective of this dissertation is to design a prac- 
tical adaptive engine control using an output regulator formulation. 

Implicit in the second objective is the need for a usable dynamic 
engine model. Such a model must be of reasonably low order while ac- 
curately predicting turbojet engine dynamics. As previously mentioned 
some work in this area has already been accomplished . However , no 
attention has been given to the important considerations of state var- 
iable selection, model order, and model structure. The usual proce- 
dure is to select a priori the order, structure, and states, identify 
a model, and verify the model. If the verification test fails another 
selection of order, structure, and states may be made and the process 
repreated until a satisfactory result is obtained (assuming a satis- 
factory result is possible using the given data and the verification 
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test). Even if the model verification is satisfactory, questions about 
the validity of the state variable selection can remain. The process 
is one of trial and error and it may be time consuming. In response to 
this need the third objective of this research is the identification of 
a low order model of turbojet engine dynamics by a technique that re- 
quires a minimum of a priori assumptions about system order and struc- 
ture. 

In summary the three objectives of this dissertation are 

(1) To develop a sampled-data (discrete) engine control algorithm 

that Incorporates adaptive sampling for efficient on-board 
computer utilization 

(2) To design a practical adaptive engine control using an output 

regulator formulation 

(3) To identify a low-order model of turbojet engine dynamics by 

a technique that requires a minimum of a priori assumptions 
about model structure and order 

To solve the problems associated with these objectives the remain- 
ing chapters are organized in the following manner. Chapter II dis- 
cusses the physics, basic control principles, and the computer simula- 
tion of a single-spool turbojet engine. The topic of Chapter III is 
the third research objective. In particular a technique by Tse and 
Weiner t (1973) is applied to the identification of turbojet engine 
dynamics. The technique requires a minimum of a priori assumptions 
and can handle stochastic output data. A model is determined for the 
turbojet engine described in Chapter II using realistically simulated 
data. In Chapter IV a digital adaptive control scheme is developed to 
Jointly satisfy the first and second objectives delineated above. 
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First, the optimal discrete output regulator problem is posed and 
solved using Lagrangian techniques for the time-invariant case. Next, 
an adaptive sampling law is developed. Finally, the optimal discrete 
output regulator and the adaptive sampling law are combined to form 
the adaptive digital control law. The adaptive digital control law is 
applied to a linearized fifth order model of a twin spool engine and 
the results simulated on a computer to evaluate its control effective- 
ness. In Chapter V the adaptive digital control scheme is applied to 
the turbojet engine simulation of Chapter II using the model developed 
in Chapter III. Results are simulated for various engine accelera- 
tions using different feedback control arrangements. Finally, this 
dissertation is concluded with a summary of results and recommenda- 
tions for future research. 



CHAPTER II 


TURBOJET ENGINES 

Turbojet engines are a conunon element In today’s modem commer- 
cial and military aircraft. Therefore, the operation and control of 
these engines is of great practical importance. This chapter dis- 
cusses the physical characteristics computer simulation, and some of 
the control concepts and requirements of turbojet engines. 

2.1 Physical Characteristics 

The purpose of a turbojet engine is to develop thrust by impart- 
ing momentum to a propellant fluid. In a turbojet this is accom- 
plished by continuously extracting, compressing, heating, and expand- 
ing air from the atmosphere. In addition to acting as the propellant 
fluid, the air also acts as the working fluid in a thermodynamic 
process . 

An ideal turbojet engine can be represented thermodynamically as 
a Brayton cycle on a classical temperature-entropy diagram (see 
Fig. 2.1). The individual processes that comprise this cycle are 

1- 2 Reversible, adiabatic (isentropic) compression between 

minimum and maximum pressures 

2- 3 Heat addition at constant maximum pressure 

3- 4 Reversible, adiabatic (isentropic) expansion between maxi- 

mum and minimum pressures 

4- 1 Heat rejection at constant minimum pressure 

In the turbojet air drawn from the atmosphere is compressed, heated. 


12 



13 



S- Entropy 


Figure 2.1. - Ideal turbojet cycle. 
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expanded, and discharged to the atmosphere by the internal engine com- 
ponents forming the continuous cycle. The internal component arrange- 
ment of a single-spool turbojet engine is shown schematically in Fig- 
ure 2.2. In flowing through these components the air undergoes sev- 
eral processes. The air is 

a-1 Brought, from far upstream of the engine, to the inlet with 
some acceleration or deceleration. Normally, this is an 
isentropic process. 

1- 2 Decreased in velocity by the inlet diffuser 

2- 3 Compressed in a dynamic mechanical compressor 

3- 4 Heated in the conbustor by mixing and burning fuel in the 

air 

4- 5 Expanded throu^ a turbine to obtain power to drive the 

compressor 

5- 6 Accelerated and exhausted throu^ the exhaust nozzle 
These processes are represented on the temperature -entropy diagram of 
Figure 2.3. In this diagram the increase in entropy due to irreversi- 
bilities are considered for each process. The effectiveness with 
which a turbojet generates thrust by these processes is highly depend- 
ent on individual component performance, the physical matching of the 
compressor and turbine, and the engine operating environment. 

The engine operating environment for a single spool turbojet 
engine (with fixed geometry) is determined by the engine fuel flow 
rate and the pressure, temperature, and velocity of the incoming air- 
stream. The airstream velocity can be given by its Mach number (ratio 
of fluid velocity to velocity of sound in that fluid) and the air- 
stream temperature and pressure are determined by the engine operating 
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Figure 2.3- “ Typical turbojet cycle. 
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altitude. Thus, a complete, but not unique, set of independent vari- 
ables that specify engine performance are Mach number, altitude, and 
fuel flow rate. Other examples of complete sets of independent vari- 
ables could be obtained by replacing engine fuel flow rate with engine 
rotational speed or the mass flow rate of air. In each case, however, 
the independent variables define an engine operating or steady-state 
point . 

Component performance characteristics are normally presented 
graphically as "maps" in terms of component pressure ratio, a rota- 
tional speed parameter, adiabatic component efficiency, and a mass 
flow rate parameter. Examples of compressor and turbine maps are 
given in Figure 2.4. The curve, denoted as ''surge line" in Figure 2.4, 
represents the boundary of stable compressor operation. Operation 
below this stability boundary is essential for satisfactory engine 
performance. 

The matching of compressor and turbine performance is a straight- 
forward problem. The turbine mass flow must equal combustor fuel flow 
and compressor airflow, and the power supplied by the turbine must 
equal that demanded by the compressor. Normally, the compressor oper- 
ates near its peak efficiency throughout its operating range when a 
desirable match is achieved. The locus of steady-state matching con- 
ditions, called an operating line. Intersects the centers of the con- 
stant compressor efficiency contours and is shown schematically in Fig- 
ure 2.5. In Figure 2.5 lines of constant temperature ratio, Tj^ (turbine 
inlet to compressor inlet} have been plotted to aid in the future dis- 
cussion of the Important problem of turbojet acceleration. A more com- 
plete discussion of jet engine turboraachinery is given by Hill and 






Compressor pressure ratio 
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Peterson (1970) . 

2.2 Computer Simulation 

Mathematical simulation techniques for turbojet steady-state and 
dynamic behavior are very useful since they enable engine dynamics and 
controls problems to be studied without endangering a valuable engine. 
Techniques for simulating an engine on analog, digital, and hybrid 
computers are fairly common and have been reported by several authors, 
e.g., Saravanamuttoo and Fawke (1970), Seldner et al. (1971), Sellers 

and Teren (1974), and Szuch (1974). 

) 

The digital simulation of a single spool turbojet used as the 
data source for the research of this dissertation was converted to a 
digital simulation from an analog simulation developed by Seldner et 
al. (1972), This simulation incorporates experimentally determined 
compressor stage data, experimentally determined lumped turbine data, 
and a real gas combustion model. The dynamics are represented by for- 
mal one-dimensional inviscld continuity, momentum, and energy approxi- 
mations to unsteady intra-compressor stage, combustor , and exhaust 
nozzle conditions. 

2.3 Control Requirements 

In general the basic control requirement of a turbojet engine is 
the determination of engine fuel flow such that (1) the engine is ac- 
celerated without violating physical engine constraints from one oper- 
ating point to another, or (2) the steady-state operating point is 
maintained in the face of external disturbances. Maintaining a 
steady-state fuel flow schedule presents no major control problems. 
Engine acceleration, on the other hand, is a more difficult control 
problem. 
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Suppose, for example, the engine is to be accelerated from oper- 
ating point A to B shown in Figure 2.5. A sudden increase in fuel 
flow causes a sudden rise in (the ratio of turbine inlet to com- 
pressor inlet temperature) before the turbomachinery has a chance to 
accelerate. Thus, the operating point moves up a constant speed line 
toward the surge line. The compressor, therefore, is moving closer 
to a region of unstable operation and possible physical damage. Addi- 
tionally, a high turbine inlet temperature may be physically damaging 
to the turbine rotor blades. Still another physical engine constraint 
is the maximum engine rotational speed. Each of these control prob- 
lems requires that fuel flow be carefully limited during accelera- 
tions . 

A basic fuel control for a single-spool turbojet is given in Fig- 
ure 2.6. This control generates a fuel flow as a function of throttle 
setting that maintains engine constraints during accelerations and 
establishes a steady-state fuel flow schedule. Several basic control 
techniques are used that are common to most fuel— speed engine con- 
trols. In this particular example, the throttle setting determines 
both a steady-state and an acceleration limit fuel flow. A compen- 
sated fuel flow error term is added to the steady-state fuel flow and 
this sum is compared to the acceleration fuel flow limit. The smaller 
of these values is compared to a lower limit (minimum fuel flow) and 
an upper limit (maximum feedback fuel flow) . If either limit is ex- 
ceeded, that limit value is selected as the engine fuel flow input. 

If neither limit is exceeded, the fuel flow selected in the first com- 
parison is the engine fuel flow input. The resultant engine speed 
output is used to generate both the fuel flow error and the maximum 



Acceleration 
fuel flow 
limit schedule 


Steady-state 
fuel flow 
schedule 

Throttle 

setting 


Figure 2.6. - Basic turbojet 
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feedback fuel flow terms. 

The acceleration fuel flow limit schedule is designed to assure 
that turbine temperature and compressor stability constraints are not 
violated. The steady-state schedule and the speed error feedback term 
assure low steady-state speed errors. The limiter is used to insure 
that the fuel flow never exceeds an overspeed fuel flow limit or falls 
below a minimum fuel flow limit. 

Certain sophistications are added to this basic control to im- 
prove engine performance. One is the addition of flight conditions, 
such as ambient temperature and flight speed as inputs to the fuel 
control. Now steady-state, feedback, and acceleration fuel flow 
schedules become multivariate functions of flight conditions and 
their respective inputs. This additional input information yields a 
better definition of engine constraints and, therefore, better engine 
performance. Another sophistication is the use of compressor dis- 
charge pressure as a control input to reduce the sensitivity of the 
control to surge producing disturbances. By normalizing or correcting 
fuel flow by a desired compressor discharge pressure, surge producing 
disturbances in the actual discharge pressure cause reductions in the 
engine fuel flow and a return to stable compressor operation. 

Figure 2.7 shows a possible engine rotor velocity trajectory and 
the resultant control modes when the basic fuel control of Figure 2.6 
is applied to a turbojet engine. The same trajectory and control 
modes are shown on a compressor map in Figure 2.8. For the trajec- 
tory the engine is initially at minimum rotor speed when a throttle 
setting change occurs that commands an engine acceleration. The dia- 
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gram shows the engine rotor acceleration, overshoot, and steady -state 
plateau. 


1 - Feedback control - minimijm steady-state speed error 

2 - Scheduled override - minimum acceleration limit 

3 - Feedback override - maximum speed limit 



Figure 2»7. - Engine speed trajectory 


1 - Feedback control - minimum steady- at ate speed error 

2 - Scheduled override - maximum acceleration limit 

3 - Feedback override - maximum speed limit 

- Feedback control - minimum steady-state speed error 



Flov rate parameter 
Figure 2,8. - Engine speed trajectory 



CHAPTER III 


ENGINE ^K)DEL IDENTIFICATION 

During the past two decades significant theoretical developments 
in the area of automatic control (optimal regulator control theory, 
Pontryagin’s ipaxlmum principle. Bellman's dynamic programming, etc.) 
have been made. The practical utilization of these theories requires 
the identification of system models from observed data to predict sys- 
tem response. Thus, much work in the identification field has been 
reported. Sage and Melsa (1971) and Eykhoff (1974) summarize most of 
these results. 

One goal of the research reported in this dissertation is to 
identify a usable dynamic model for a gas turbine engine. This model 
will be used in the application of modem control theory to turbine 
engine control. Work in identifying gas turbine engine models has 
been done by Otto and Taylor (1951) , Crooks and Willshire (1956) , 
Ahlbeck (1966), Chen (1972), Michael and Farrar (1973), and Mueller 
(1971) . These techniques assume a model order and linearize the non- 
linear dynamics about an operating point to obtain linear operating 
point models. The parameters in the models of Ahlbeck, Mueller, and 
Chen were determined from transfer ftmction analysis, Michael and 
Farrar used a least squares curve fitting technique to find the system 
parameter values. These examples show that, for control purposes, 
linearized operating point models can adequately represent the non- 
linear response of turbine engines. 
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In the previously mentioned papers on engine identification the 
model order and the model structure were assumed. Then the parameters 
were determined from available system data for this assumed structure. 
If these initial model assumptions are incorrect, the result may be a 
model that is too complex or too simplistic for design and control 
purposes. One of the contributions of this dissertation is to apply 
an identification technique that will not only identify the param- 
eters of a suitable linear operating point model of a turbojet engine, 
but also determine the appropriate model order from accessible engine 
data. Such a technique eliminates the need for initial model order 
and structure assumptions and has been developed by Tse and Welnert 
(1973) for autonomous linear discrete systems. 

The remainder of this chapter describes the Tse and Weinert iden- 
tification technique and its application to the identification of 
autonomous linear operating point models of the single-spool turbojet 
described in Chapter III. To complete the model, parameters relating 
system control variables to engine response are identified by gradient 
search. Finally, the engine response as predicted by the completed 
model is compared with the actual response of the simulation. 

3.1 Tse and Weinert Identification Technique 

Tse and Weinert have developed a method that identifies from out- 
put data a constant, multivariable, stochastic, linear system which 
has xmknown dimension, system, matrices, and noise covariances. A 
general stochastic model is not identifiable from steady-state output 
data since the output data determine an equivalence class of systems 
(Tse and Weinert, 1973) . The systems in this equivalence class have 
steady-state Kalman filters with the same impulse response and inno- 
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vations covariance (Tse and Anton» 1972; Tse and Welnert, 1973). If 
the system matrices are chosen in a certain unique canonical form 
(Weinert and Anton, 1972), it is possible to obtain consistent esti- 
mates of the Kalman filter parameters. The procedure for consist- 
ently estimating the Kalman filter parameters, transition, observa- 
tion, gain, and noise covariance matrices, and noniteratively estimat- 
ing system order is described below. 

Problem statement . Consider the linear discrete system repre- 
sented by 

x(k + 1) " Ax(k) + w(k) 

(3.1.1) 

y(k) = Cx(k) + v(k) 

where xeR’^, yeR“, and w(k) and v(k) are zero-mean Gaussian noises 
with coArairiances 

E{w(k)w’(j)} = 

E{v(k)v'(j)} = (3.1.2) 

E{w(k)v'(j)) = 

The unknown vector of parameters is 

0 = {xq,A,C,W,V,D} (3s1.3) 

The object is to identify Q using the observed output data 
= {y(l),y(2), . . .,y(N)}. 

Tse and Anton (1972) have shown that using the observation data 
the appropriate parameter vector to identify is ■ {A,C,B,Q} 
where the steady-state Kalman filter associated with (3sl*l) is 
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x(k + 1/k) = Ax(k/k - 1) + Bw(k) 
u(k) » y(k) - Cx(k/k - 1) 

and 

E{y(k)u'(j)} = 


(3.1.3) 


(3.1.4) 


The variable u(k) is the innovation process (Kailath, 1970) and B 
is the steady-state Kalman filter gain given by 


B = (APC’ + D)Q 


(3.1.5) 


P = APA' + W - BQB’ 


(3,1.6) 


Q » CPC' + W 


(3.1.7) 


P is the steady-state one-step prediction error covariance. The 
assumptions on 6^ are 

(1) A is stable, 

(2) (A,C) is an observable pair. 

(3) (A,B) is a controllable pair. 

(4) (A - BC) is stable. 

(5) The system dimension n is finite and unknown. 

(6) The effect of any initial condition on the system has died 
out (i.e., the system is in steady-state). 

System canonical form . A canonical structure for (A,C) derived in 
Weinert and Anton (1972) is summarized as follows. Let c| be the i^^ 
row of C. Consider the rows of the observability matrix in this 


order : 
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O* 

f > 

* * * * 


(3.1.8) 


Pi 

Let p^ be the smallest non-negative integer such that is 

linearly dependent on the vectors in all the preceding rows of (3.1.8) 
Then 


m 


i=l 


(3.1.9) 


By the definition of {p . , there is a unique set for 

X X“X 

i - 1,2, . • . ,m such that 


cl*^' 


i)D 

j«l k=0 


if p. > 0 
1 


(3.1.10) 


i-1 Pi“l 

“ Zy Zy 

j=l k-0 


if p^ = 0 


(3.1.11) 


Given the set ^Pi»^ljlt^^» ® canonical form for (A,C) is uniquely speci- 
fied by 


A « 


^11 


ijml' * * 


(3.1.12) 
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r 1 

0 1 I 


A — 

1 


A ■ • “ 

11 



(pi p^) 

fiiO’^iil’ • 

• '*®ii,pi-l^ 

A 

0 





(Pi Pj) 

i > j 

_fijO»^ijl’ ’ 



and 


~ [ 0 , • • • • • » 0 ] 


Pi 


> 0 


(3.1.X3) 


(3.1.14) 


(3.1.15) 


where the 1 in is in column 1 + + P 2 + . . . + Pi_i> 

if Pj^ = 0 

*^i “ ^^ilO» • • * »^il,Pj^-l*^i20 ®i2,p2-l» * * *» 

^i,i-l,p^_j^-l’° (3.1.16) 


The matrix B has no special form but its entries and 
uniquely determined by (Weinert, 1973). Therefore to identify 
the system dimension n, the system transition matrix A, and the ob- 
servation matrix C, it is sufficient to estimate Tse and 

Weinert estimate these parameters using time series analysis. 

Time series analysis . As stated above the estimation of 
{p^j^ljj^} is sufficient for the identification of n. A, and C. Time 
series analysis is used to accomplish this estimation. Consider ti.^ 
system 

z(k + 1) - Az(k) + By(k) 

(3.1.17) 

y(k) » Cz(k) + y(k) 
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which is equivalent to (3.1.3) when 

E{p(k)y*(j)} = (3.1.18) 

Now, let 

R(a) * E{y(k + a)y*(k)} a = 0,1,2, . . . (3.1.19) 

In Appendix A it is shown that the set ca^ determined 

from the elements of R(a) . This is accomplished as follows . Let 
r^j(c) be the i,j^^ element of R(a) . 

R(a) = lr^j(o)l (3.1.20) 

Define 


k 

\ = ( 3 . 1 . 21 ) 

i-1 

Note that = n, where n is the system dimension. Also define 
the following vectors 

+ 2 ), . . ( 3 . 1 . 22 ) 


and 

8^ = (zi j Z2 j . . .1 

^ii,Pjj^-i^ 

Also, define the ideutlf lability matrix as 


(3.1.23) 

(3.1.24) 
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'2I 


I ,j ^i+1, j 

I 


'z.l 

I ^ 

f 




(3.1.25) 


where 


«J(W 


j (1) »r (2) j , . ,,rj^j(k) 
ri^(2),rij(3), . . .,rj^j(k + l) 


r^jCk) ,r]^j(k + 1), . • . ,r (2k - 1) 


l'~'1^2f . . .y m end 


(3.1.26) 


+ 1) > • • * ^ 

r^j(L± + 2). . . .,r^j(L^ + ?£ + ^> 


(L^ + k) , « • + p^ + k — 1) 


(3.1.27) 


The index j may take any integer value between. 1 and m in equa- 
tions (3.1.22) to (3*1.27). 

Tse and Weiner t (1973) show that (see Appendix A) 
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(3.1.28) 

Defining 

= Determinant ^(k) (3.1.29) 

The values for (p^} are found by testing d^(k) , k = 1,2, . . . until 
<i|(q) = 0, in which case * q - 1. The parameter 8^ is then found 
by solving (3.1.28). 

Since only is available, R(a) will not be known exactly. It 

can be estimated, however, as follows: 

n-g 

R(g) « ^ y~^y(k + g)y*(k) (3.1.30) 

k-1 


Using estimated values in (3.1.22) to (3.1.29) p^ is obtained by se- 
lecting some threshold e and testing d^(k) , k “ 1,2, . . . until 
|d^()l)| < e. Now p^ « ^ - 1. The parameter 8^ is the solution to 


^ A 


n - 


(3.1.31) 


From (3.1.19), (3.1.30) and Parzen (1967), it can be shown that 
r^^(g) is an asymptotically (N ->• <») unbiased, normal, and consistent 
estimate of r^j(a). Results in Mehra (1971) can be used to show 
that 8^ is an asymptotically unbiased, normal, and consistent esti- 
mate of 8^, given that p^^ = p^. 

3.2 Application to Engine Simulation 

This section discusses the development of linearized operating 


point models of engine dynamics, and the application of Tse and 
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Weinert’s method to the identification, of these operating point models. 
The digital simulation of the turbojet engine described in Chapter III 
is the data source for this dissertation. Linear operating point 
models are determined from this simulation using only inputs and out- 
puts realizable at an acttial engine. 

Linearized operating points . It is assumed that the dynamical 
input-output-state relations for the ordinance engine are given by the 
nonlinear, vector differential and algebraic equations 


* 

X = f(x,u) + C 
y * y(x,u) + Y 


( 3 . 2 , 1 ) 


The vector xeR*^ represents the state of the system, the vector ueR^ 
is the system input, and the vector yeR™ is the system output. The 
vectors 5 and y are Gaussian white noise vectors with unknown sta- 
tistics. The functions f(. , ,) and g(. , .) are assumed continuous 
and twice differentiable in their arguments. The vectors x, u, 
and Y functions of time. 

Expanding the functions f(. , .) and g(. , .) in a Taylor 

series about the steady-state operating point (x ,u ) results in the 

ss s s 

following system equations vdiere f(x„_,u _) ® 0. 

s s ss 


X +6 
ss X 


f(x u ) + 


fix 


x«x, 

u=u, 


ss 

ss 


^ 3f(x>u) 


du 


fill + HeO*T. + ^ 




( 3 . 2 . 2 ) 
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y 

-'ss 




5x 


x=x 

u=u 


ss 

ss 


^g(x>u) 

3u 


3C==Xss 

u=Uss 


6u + H . O.T , 


+ Y 


(3.2.3) 


Dropping the higher order terms (H.O.T.) and defining 


F = 


3£(x,u) 

dx 


x=Xss 

«=«ss 


G = 


3f (x,u) 
3u 


x=x, 

u=u, 


ss 

ss 


(3.2.4) 


H 


3g(x,u) 

3x 


x=*ss 

u=uas 


E 


3g(x,u) 

3u 


x=x 

u=u 


ss 

ss 


and simplifying equations (3.2.3) gives the first order or linearized 
approximations to equation (3.2.1) 


6x = F6x + G6u + ? 
6y = H6x + E6u + Y 


(3.2.5) 


Mathematically (3.2.5) constitutes a set of linear, constant 
coefficient, multivariable, stochastic differential and algebraic 
equations . If this linearization procedure is accomplished over a 
sufficiently large number of steady-state operating points, engine 
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dynamics can be approximated over the entire operating range of the 
engine. 

For the single-spool turbojet described in Chapter III, an oper- 
ating point is uniquely specified by engine rotor speed, flight Mach 
number, and flight altitude. In this dissertation the engine is 
assumed to operate at a sea level static condition. This standard 
test condition specifies the Mach number and the flight altitude. 
Therefore, the engine operating point, and consequently the operating 
point models, vary only with engine rotor speed. 

The linear, time-invariant, operating point models identified in 
this chapter are used to construct a composite engine model and to 
generate output feedback gains. The composite engine model will be 
used to verify the identification results of this chapter by compar- 
ing composite model and simulation dynamics. The feedback gains will 
be generated from individual operating point models and combined into 
the adaptive digital control described in a subsequent chapter. 

Application of Tse and Weinert^s method . The model equations , as 
determined in the previous section, for an operating point are 

6x F6x G6u + C 

(3.2.5) 

6y = H6x + E6u + Y 

Because Tse and Weinert^s method requires a discrete model, a discrete 
version of equation (3.2.5) is required. Following the procedure of 
Appendix D a discrete version of (3.2.5) is 
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6x(k + 1) = A6x(k) + B6u(k) + C^(k) 

( 3 . 2 . 6 ) 

6y(k) = C6x(k) + D5u(k) + Yj^(k) 

if 6u(t) is assumed constant over the sampling period, T, and 


* FT 

A = e 

(3.2.7) 

OT 

B = / dt 

(3.2.8) 

f 

C = H 

(3.2.9) 

D = E 

(3.2.10) 


Also, it is assumed that Cj^(k) and Yj^(k) are discrete Gaussian white 
noise vectors with zero mean and unknown covariances. 

E(c^(k)}=0 

( 3 . 2 . 11 ) 

E{Yi(k)} = 0 

The digital turbojet simulation used in this dissertation does 
not Include possible engine noise sources such as random variations in 
the compressor inlet conditions or the combustion process. To further 
simulate a real engine and to facilitate identification, zero mean 
white noise was introduced into the simulation. This was accomplished 
simply by adding a Gaussian random number to the value of rotor speed 
at each iteration of the simulation. 

Eight outputs (ra = 8) were selected for the identification pro- 
cedure , They are 


y* ~ ^ ^ *^C *"^T *^T 


( 3 . 2 . 12 ) 
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where the variables are defined as follows 
S rotor speed 

compressor discharge temperature 
Pg compressor discharge pressure 
T^r nozzle inlet temperature 
nozzle inlet pressure 
turbine inlet temperature 
Pj turbine inlet pressure 
engine thrust 

The first five variables are all readily measurable. 

Initially a sample transient was simulated to determine approxi- 
mate dynamics. Since the Tse and Weinert procedure requires data that 
is representative of the system as t some estimate of system 

time constants is needed to insure that the interval of data taking is 
sufficiently long. For a step change in engine fuel flow, w^, the re- 
spective change in rotor speed, S, is shown in Figure 3.1. The rotor 
speed time constant (typically the engine's largest time constant) is 
approximately 1 second. Thus data taken over a 10 to 15 second inter- 
val will adequately approximate data taken over an infinite interval. 
The maximum frequency of Interest is assumed to be 5 hertz. Thus, the 
data sampling period is T = 0.1 second as determined by the sampling 
theorem (Shannon, 1949). 

The Identification algorithm requires the vector <5y(k). D._I-wn- 
ing y(k) as 

y(k) « + 6y(k) ( 3 . 2 . 13 ) 

NOM 


where v 


is some constant nominal vector and y(k) is available 
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Time (sec) 

Figure 3.1. - Response of rotor speed to a step 
change in fuel flow. 
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from the simulation. Also, the Identification algorithm does not 
allow for a deterministic input. Thus the input vector u(t) will be 
held constant which Implies that 

6u(k) = 0 (3,2,14) 

Then 

6x(k + 1) = Ax(k) + Cj^(k) 

(3,2.15) 

6y(k) = Cx(k) + Yi(k) 

Since steady-state is assumed, the effect of x(0) has been eliminated. 
Then since 


6x(k) 


6x(0) + \(i) 

i=l 


and 


E{6x(k)} = E< ^A^"^Cj^(i)> = 0 

.i=l 


then 


E{6y(k)} = E{C6x(k) + Yi(k)> = 0 


From (3,2,18) and (3,2.13) 


E{y(k)) - 


(3.2.16) 


(3.2.17) 


(3.2.18) 


(3.2.19) 


Since the data interval is finite, the vector Fj;|qjjj will be esti- 


mated as 
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N 

y *l\^y(l) (3.2.20) 

i=l 

where N is the number of data points. Then 

N 

6y(k) = y(k) - I ^y(i) 
i=l 

The time histories {6y(k) are the data required by the Tse and 
Weiner t algorithm. These time histories were simulated at four dif- 
ferent operating points. Figure 3 . 2 shows typical trajectories for two 
components of the output vector defined in equation (3.2.12) when the 
simulation is disturbed by noise. 

The first step in the algorithm is the determination of the model 
structure. The model structure is defined by the parameters 
The parameters are estimated as outlined previously by determining 
when the determinant of the identification matrix (3.1,25) falls be- 
low a certain threshold value or exhibits a sharp decline in value. 

From the definition of the observation matrix C it can be seen 
that at least one output variable can be made a state variable of the 
identified system. Since much of the control work to follow will be 
centered around rotor speed, and since the operating point is a func- 
tion of rotor speed, rotor speed was the output variable initially 
selected as that state variable. 

To determine the structure parameter pj^ the correlation matrix 

A 

estimate [R(o) ] was calculated from the observed data for the four 
chosen operating points. Then the Identification matrix was 

formed from the elements of R(o) . The determinant of 4>^(k) , i = 1 



S-Hotor speed (lOOO rev/min) 



Time (sec) 

(a) Rotor speed. 


Figure 3«2. - Typical output trajectories for a 
Gaussian disturbance. 



Time (sec) 

(b) Turbine inlet temperature. 

Figure 3.2. - Concluded. Typical output trajec- 
tories for a Gaussian disturbance. 
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for various values of k and j is given in Figure 3.3, Typical 
plots for other j values give similar results* From the plots and 
the criteria mentioned above it appears that a first order system 
would adequately model the dynamics of the rotor speed output vari- 
able* However, a second order model was chosen to give greater con- 
trol design flexibility and more accurate prediction of system re- 


sponse* 

Next the matrices *^(k) , i = 2, . . .,8 were calculated for 
various values of j and k. For each of these matrices the respec- 

A 

tive estimate of the structure parameter, p^, was equal to zero. This 
indicates that each corresponding output is a linear combination of 


the state variables. 

The structure of the model now becomes using equations (3.1.12) 
to (3.1.14) 



since pj^ = 2. From equation (3.1.16) 

r- ~l 


1 

0 

^210 

®211 

^310 

^311 

^410 

®411 

^510 

^511 

^610 

^611 

®710 

®711 

^SIO 

^811 


( 3 . 2 . 22 ) 


(3.2.23) 
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since « 0, i »= 2,3, . . . ,m; m = 8. Once the model has been 
parameterized, the next step Is to estimate the values of these 
parameters. These estimated parameter values are determined by equa- 
tion (3.1.31) at each operating point. 

Four operating points were selected to approximate the engine’s 
operating range. Designating the design speed of 36 960 rev/sec as 
100 percent rotor speed, the four operating points selected corre- 
spond to 80, 90, 100, and 104.5 percent rotor speed. This set of op- 
erating points represents the endpoints (80 and 104.5%), the design 
point (100%), and an intermediate point (90%) of the operating range. 
The identified discrete system matrices, A and C, for each operat- 
ing point are given in Figure 3.4. Since the continuous system 
matrices will be required in a subsequent chapter, F and H are 
calculated from their discrete counterparts, A and C, using the 
concept of a logarithm of a matrix (Gantmacher, 1959) and listed in 
Figure 3.5. Appendix B shows the procedure followed to transform 
the discrete system to a continuous one. 

3.3 The Control Matrices 

The discrete model for the system is given as 


6x(k + 1) = A6x(k) + B6u(k) + ^j^(k) 
6y(k) = C6x(k) + D6u(k) + Yj^(k) 


(3.2.6) 


Tse and Weinert’s algorithm was used to determine A and C in a 
specified canonical form. For this canonical form the control matrix 
B and the direct link matrix D will have no special form. Thus B 
has n X q elements and d has ro x q elements that must be identi- 


fied 



^ Speed Matrix 


Matrix elements 



Ln 


Figure 3.4. - Identified values of A and C vs. ^ speed-discrete model. 



4^ 


^ speed Matrix Matrix elements 


80 

F 

-15.4258 16.6030 

-5.8816 5.0482 

H’ 

, 

Same as C 

; in Figure 3 A 

90 

F 

-15.4915 17.0706 

-5.8021 4.7002 

h' 

Same as C' 
in Figure 3.4 

100 

F 

-17.0582 19.4081 

-5.0076 3.5108 

H’ 

Same as C’ 
in Figure 3.4 

104.5 

F 

-15.6751 17.7771 

-5.6552 4.2211 

H 

Same as C 
in Figure 3.4 


Figure 3.5, - Identified values of F 
and H vs, ^ speed- continuous model. 
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The order of the system has been determined, n = 2. The control 
variable of the simulation is fuel flow, thus q = 1 and the B 
matrix has two unknown elements . Since the number of output variables 
is 8, the number of unknown elements in the D matrix is 8. 

These ten elements were determined at each operating point by a 
simple gradient search procedure. First a control input was selected 
as 

6u(k) = <5w^(k) = A[sin(0.1 wkT + 

+ sin(wkT + + sin(10 wkT + (3.3.1) 

where A = 0.00667, w = 2, = 0.1, - 0, and = 1. This con- 

trol input was applied to the nonlinear simulation and y(t), the out- 
put vector, was calculated for a 5-second interval. Then an Initial 
guess, and D^, for the matrices B and D was chosen and the 
equations of (3.2.6) simulated to give y^^(k) using B = and 
D = D^. The trajectory y(t) was compared at the appropriate sampling 
points to the output estimate 

y(k) * yjj + 'Syi(k) (3.3.2) 


The squares of the errors were summed to form a cost function J^. 

N 

*7^ = ^ [y(k) - yi(k)]'[y(k) -y^(k)J (3.3.3) 

k=l 


Another guess for B and D, called 
by perturbing one element of the 
^ERR' estimated trajectory 


®i+l ^i+1’ determined 

or matrix by a small amount 

(k) again was simulated and com- 
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pared to y(t) . The cost function ^£^2. calculated according to 

(3.3.2). The gradient of the cost function with respect to one ele- 
ment of the matrix B or D is estimated by 




(3.3.4) 


where is one of the elements of the B or D matrix. A new 

choice for b ^^^2 then be found according to 

8J. 

'“i+2 ■ >>1 + 


This procedure is followed successively for the first and second ele- 
ments of B and the eight elements of D until the gradient for 
each element becomes small. The values for B and D determined by 
this gradient search procedure are listed in Figure 3.6. 

The continuous counterpart, G, to the B matrix can be found 

from 


B 



(3.3.6) 


by integrating and solving for G. 

G = (A - I)“^FB (3.3.7) 

The G matrices for the four operating points are also listed in 
Figure 3.6. The continuous matrix E is equal to the discrete matrix 
D and is listed in Figure 3.7. 



Figvire 3.6. - Identified values of 
B and G vs. ^ speed. 







^ Speed Matrix 


Matrix elements 
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E-' 

-12227.78 -15.25 12.47 3508.81 2.598 3991.66 13.96 -11703.98 

90 

E' 

-10573.77 -33.38 4.550 2920.59 2.970 3427.02 9.0856 114.41 

1 

100 

! 

E' 

-18929.90 -223.59 -26.322 2535.71 -.9140 2870.89 -23.632 -117.678 

10J4.5 

E* 

-14468.60 -203.73 -19.96 2299.11 .5885 2632.15 -20.069 -43.53 


Figure 3.7. - Identified values of D and E vs. ^ speed. 





51 


3.4 Model Verification 

From the operating point models obtained at four engine rotor 
speeds, a continuous function of speed was calculated by linear inter- 
polation for each model parameter . The model equations become 

5x(t) = F(S)6x(t) + G(S)6w^(t) (3.4.1) 

6y(t) = H(S)6x(t) + E(S)6w^(t) 

where the vector y(t) is as defined in equation (3.1.12) and S is 
rotor speed. 

The composite model was simulated on a digital computer. A 
block diagram of the composite model is given in Figure 3.8. Note 
that the nominal fuel flow, is found by time averaging the 

engine fuel flow, w^. Since steady -state fuel flow determines an 
operating point, the nominal output vector, defined 

as the steady-state engine output if a constant fuel flow, is 

supplied to the engine. Thus, the nominal output vector can be 
scheduled as a function of nominal fuel flow. 

A test input was selected as a combination of three basic con- 
trol inputs, the step, ramp, and parabola. This test input, shown in 
Figure 3.9 was applied to the engine simulation and to the composite 
model. The engine simulation output, the composite model output, and 
their difference (error) are plotted in Figure 3.10. Also a listing 
of the average error and error variance for each output is given in 
Figure 3.11. The normalized average error for each output was less 
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("b) Compressor discharge temperature* 


Figure 3- 10* - Simulation and composite model trajectories. 



Temperature (^R) 
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Compressor discharge temperature error. 


Figure 3»10. ^ Simulation and composite model trajectories. 


-Pressure (psi) 
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(c) Compressor discharge pressure. 

Figure 3-10* — SimuXation and composi'te model “trajectories • 
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(cc) Compressor di 


Figure 3*10. - Simulation 
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Output 

Average 

error 

Variance 

Normalizing 

value 

Normalized 

average 

error 

Normalized 

variance 

S 

260.960 

138860 

35000 

.0077 

.0106 

Tc 

4.396 

14.06 

810 

.0054 

.00463 

Pc 

.8092 

.3057 

51..53 

.0157 

.0107 

Tz 

10.47 

2484 

1283 

.0082 

.0389 

Pz 

.1279 

.01484 

20 

.0064 

.0061 

Tt 

15.47 

3345 

1544 

.0100 

.0375 

Ft 

.7351 

.2497 

46,51 

.0158 

.0107 

Pz 

7.615 

86.05 

327,19 

.0233 

.0284 


Normal average error = 


Average error 
Normal value 


Normal variance = 


Variance 


(Normal value) 


2 


Figure 3.11. - Average error values and variances for output 
error trajectories. 
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than 3 percent. This error is considered insignificant for control 
purposes and thus the model passes the verification test. Fig- 
ures 3.12 and 3.13 are plots of as determined by the 

composite model. 

It is remarked that a more accurate model, if required, could be 
attained by three methods. First, the estimates obtained by the Tse 
and Weinert method could be enhanced by using more data in the identi- 
fication procedure. More data implies either longer time history 
intervals or several shorter time history intervals with the result- 
ing parameter estimates averaged together. Second, a more accurate 
identification technique that requires good a priori information could 
be applied. Using the structure and the parameters determined by the 
Tse and Weinert method as the assumed model structure and initial 
parameter conditions, the more accurate identification method would 
have adequate a priori information for a fast and accurate conver- 
gence of parameter estimates . Such methods include quasillnearlza- 
tlon, invariant imbedding, and sequential identification (Sage and 
Melsa, 1971) . Third, the data sampling period could be reduced to 
include higher frequency components in the identification data. Con- 
sequently, the model would be accurate over a larger frequency range 
and reduce the error caused by high frequency elements. 





CHAPTER IV 


ADAPTIVE DIGITAL CONTROL 

This chapter discusses the development of an adaptive digital con- 
trol scheme to be used for jet engine control. This adaptive digital 
control scheme must jointly satisfy the first and second research ob- 
jectives stated in Chapter I. These objectives are to (1) develop a 
engine control algorithm that xncorporates adaptive sampling 
and (2) design a practical adaptive engine control using an output 
regulator formulation. 

In a digital control system there exists a sample and hold 
operation that quantizes the continuous analog signal of the continu- 
ous controlled process into a discrete digital signal capable of being 
processed by a digital computer. One of the major engineering consid- 
the design of digital control systems is the sampling 
frequency. If this frequency is too small, system instabilities will 
occur. If the sampling frequency is too large the computer will be 
required to process more xnformation than xs really requxred for 
adequate control. An alternative to choosing a fixed sampling fre- 
quency for the digital controller is adaptive sampling. Adaptive 
sampling increases sampling efficiency by varying the sampling fre- 
quency as a function of system parameter. Sampling efficiency is 
defined as the ratio of some quantitative measure of a system perform- 
ance to the number of sampling Instants required to achieve that 
performance. 


67 



68 


Adaptive sampling was first developed by Dorf, et al. (1962). 
Later other authors such as Gupta (1963) , Tottiovic and Bekey (1966) , 
Mitchell and McDaniel (1969), and Bekey and Tomovic (1966) developed 
alternate adaptive sampling schemes. Hsia (1974) has shown that the 
design of adaptive sampling laws could be unified under one analytical 
approach. 

Since adaptive sampling increases sampling efficiency, the same 
level of system performance can be achieved with adaptive sampling as 
with fixed frequency sampling but with fewer sampling instants. Also, 
computer usage for control purposes is directly proportional to the 
number of sampling instants. Thus adaptive sampling applied to jet 
engine control can make efficient use of on-board computers, one of 
the objectives of this research. 

The second objective of the research reported in this chapter is 
the design of a jet engine control that is simple, practical, and 
maintains good engine response. A simple, practical, and effective 
way of regulating the outputs of many systems is through constant, 
proportional feedback. To Insure good system response, constant pro- 
portional feedback laws can be designed via optimal regulator theory, 
(Kwakernaak and Sivan, 1972). Michael and Farrar (1973) have shown 
that the continuous optimal state regulator theory can be applied to 
the design of controls for a jet engine. One of the assumptions 
necessary for the work of Michael and Farrar was the availability of 
all state variables. This is not always practical in jet engine con- 
trol. Another solution is to reconstruct the states from the engine 
outputs using a Kalman filter (Kalman and Bucy, 1961) or a Luenberger 
observer (Luenberger, 1966). However, this would probably result in 
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an overly complex control. An alternate approach is to use a time- 
invariant optimal output regulator. Such a control would be both simple 
and practical since the feedback gains are constant and only output 
measurements are required. Levine and Athans (1970) state that the out~ 
put regulator will perform well for many well-behaved systems. To fa- 
cilitate the use of an on-board digital controller the optimal output 
regulator solution is required in discrete form. Thus the second ob- 
jective will be met by an optimal discrete output regulator. 

This chapter is divided Into three sections. In the first sec- 
tion, the discrete output regulator problem for time- invariant linear 
systems is stated and necessary conditions for its solution are de- 
rived. Also, an algorithm for computer solution of the necessary con- 
ditions is given. The second section describes adaptive sampling and 
the derivation of adaptive sampling control schemes. Also computer 
simulation results of the application of this control scheme to a 
linearized fifth order jet engine model are presented. 

4.1 The Optimal Discrete Output Regulator 

The optimal linear state regulator is a well known and well 
studied problem. The fundamental results are by Kalman (1960) . Sev- 
eral texts that extend the basic results in both discrete and contin- 
uous time formulations are also available, see for example that of 
Anderson and Moore (1971) and Kwakernaak and Sivan (1972). 

The optimal linear state regulator requires the full state of 
the system to determine the feedback control. Often the order of the 
output vector of practical systems is less than the order of the sys- 
tem state vector. Thus to apply the optimal linear state regulator, 
either a Kalman filter (Kalman and Bucy, 1961) or a Luenberger 
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observer (Luenberger, 1966) is often used to generate an estitnate of 
the state vector. 

An alternate approach is to design a regulator that uses only 
available outputs. Results on this specific problem have been obtained 
for systems with a scalar control by Rekasius (1967) and for multivar- 
iable systems by Levine, et al. (1970 and 1971). Mendel (1974) has 
also obtained similar results. 

For the discrete time case Mullis (1973) has developed weak suf- 
ficient conditions for the existence of a finite sequence of output 
feedback gains for which every initial state can be driven to the ori- 
gin. Ermer and VandeLinde (1972) have developed the discrete output 
regulator for the time-varying case using dynamic programming. 

In this section a sampled-data constraint is imposed on a contin- 
uous system with a quadratic cost function. The resultant discrete 
or sampled-data system will have cross weighting between the state 
and the control in the cost function. For full state feedback a 
quaidratic cost function with cross weighting can be converted easily 
to one without cross weighting by a transformation Involving the state 
and control variables. In the output feedback formulation such a 
transformation is not physically realizable since the full state vec- 
tor is not available. This section, therefore, develops the discrete 
output regulator equations when the process dynamics are linear and 
time- invariant for a cost function with cross weighting between the 
state and control vectors. An iterative algorithm is developed that 
solves the equations of the time- invariant case, and an example prob- 
lem is solved using this algorithm. 
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By analogy to the time- invariant state feedback regulator problem, 
one might expect that the optimal output feedback matrix would become 
time-invariant as the interval of control becomes semi-infinite. How- 
ever, for many well behaved systems this is not the case. Brockett 
and Lee (1967) have cited an example of a second order system with one 
output that is both observable and controllable. The system, however, 
cannot be stabilized by a constant feedback gain. Yet the system can 
be stabilized by a time-varying gain. To avoid this difficulty, 
the formulation to be described constrains the feedback matrix to be 
a constant and assumes that a constant stabilizing output feedback ma- 
trix does exist. 

Problem Formulation . Given the time -Invariant linear system 

x(k + 1) = Ax(k) + Bu(k) (4.1.1) 

y(k) = Cx(k) 

where x e y e u e Consider also the quadratic cost func- 

tion 

OQ 

J = }^x'(i)Qx(l) + u'(k)Ru(i) + 2x'(i)Mu(i) (4.1.2) 

i=0 

where it Is assumed that 


A 

is 

an 

n X n 

real constant matrix 

B 

is 

an 

n X m 

real constant matrix 

C 

is 

a 

q X n 

real constant matrix 

Q 

is 

an 

n X n 

symmetric positive semi-definate real con 


stant matrix 


M is an n X m real constant matrix 
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R is an m X m symmetric positive definate real constant 

matrix 

Q - MR“%' is positive semi-definate 
Now introducing the constraint 

u(k) = -Sy(k) = -SCx(k) (4.1.3) 

the cost function becomes 

•^“7 x' (k) IQ + C'S'RSC - MSC - C'S'M']x(k) (4.1.4) 

^ k=0 

Using a theorem from Kwakernaak and Sivan (1972) , the constrained dy- 
namic optimization problem can be converted to a constrained static 
problem. 

Theorem 1 (Kwakernaak and Sivan, 1972) . Let x(k) be the solution 
of 


x(k + 1) = Tx(k) (4.1.5) 

x(0) = Xq 

If T and the symmetric positive semi-definate matrix V are constant 
and the moduli of the characteristic values of T are strictly less 
than one, then 




Tr PEIxqX^ 


(4.1.6) 


where 

P = T'PT + V (4.1.7) 

Taking the expected value and applying Theorem 1, the cost function be- 


comes 
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J = E(J) = y THPXq} 


(4.1.8) 


with 


P = 


A(JPAq + V 


(4.1.9) 


and 


V “ Q + C'S’RSC - MSC - C'S'M' 


Aq = A - BSC 


Xq = e(xox5) 

The optimization problem now is to find a matrix S* that minimizes 
the performance index of equation (4.1.8) while satisfying the con- 
straint of equation (4.1.9). 

The problem is solved by defining the Lagrangian 


(4.1.10) 

(4.1.11) 

(4.1.12) 


X = i TrjpXo + (^-P + AqPAo + v)l'^ 


(4.1.13) 


where L is a matrix of Lagrange multipliers. The necessary condi- 
tions for optimization are 





, a)f 

= 0. 

* 

ap| 

- 0, 

1 * 



as 


Using the following formulas 


-f- {Tr[AXB]> = A’B' 
dx 


= 0 


(4.1.14) 


ax 


{Tr[AXBX ]} = A’XB' + AXB 


(4.1.15) 

(4.1.16) 


The necessary conditions become 

S* = (R + B’P*B)“^(B’P*A + M')L*C'(CL*C')”^ 
L* - aJl*A^* + Xq 


(4.1.17) 

(4.1.18) 



P* = A^’PAq + V* 
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(4.1.19) 


where 

Aq = A - BS*C (4.1.20) 

and 

V* = Q + C'S’*RS*C-MS*C - (MS*C)' (4.1.21) 

These equations represent necessary conditions for the optimal 
discrete output regulator. Both equations (4.1.18) and (4.1.19) are 
in the form of discrete Lyapunov equations. However, if S* were 
eliminated in equations (4.1.18) and (4.1.19) via (4.1.17), then equa- 
tions (4.1.18) and (4.1.19) would become a pair of coupled discrete 
Ricattl equations. Also, note that if C is square and invertible 
then these necessary conditions reduce to the steady-state equations 
for the discrete optimal state regulator. 

Computer Algorithm . A computer program was written that solves 
equations (4.1.17) to (4.1.21). The algorithm that was programmed is 
as follows. Given the equations 

(1) Si = ^R + B’PiB^'^^B'PiA + M'jLiC' [cLiC'J"^ 

(2) Ai?i+iAi - Pi+i = -Q - C'SiRSiC +MS^C+ (MSiO ' 

(3) AiLiAi - Li - -I 

(4) Ai = A - BSiC 
where now 

Xo = I 

(a) Set i = 1, Pi = 0 

(b) Select an initial stabilizing gain matrix. Si 

(c) Calculate Pi+j using (2) and (4) and Si 

(d) Stop the iterative procedure if j |Pi - Pi+il I e 
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(e) Set i = i + 1 

(f ) Calculate simultaneously and using (1) , (3 ) , (4 ) , 

and 

(g) Return to step (c) 

The subroutines used in this algorithm are entitled CLSDLP, MULT, 
DISLYP, DITORF, and RICATT. They were written in Fortran IV for the 
IBM 7094 digital computer. A listing for each subroutine is given in 
Appendix C. 

The subroutine MULT determines the product of two compatible ma- 
trices. The subroutine CLSDLP solves for the closed loop system ma- 
trix and its transpose given A, B, S, and C. Subroutine DISLYP 
solves the matrix equation 

YPV _ p = (4.1.22) 

by successive approximation given a stable matrix V and a positive 
definate matrix W, The subroutine DITORF solves equation (4.1.17). 

The subroutine RICATT is the controlling subroutine that accomplishes 
the iterative portion of the algorithm. Convergence criteria used in 
RICATT and DISLYP require a matrix norm. If Aj is an n x n matrix 
with elements aj^^, then the norm of this matrix as used in RICATT and 
DISLYP is defined, as 

lUdll = fif: layl W-l-23) 

i=l j-1 

It is remarked that an initial gain matrix, for the full order 

case (n = m) , can be obtained from an algorithm by Kleinman (1974) . 

For the reduced order (n ^ m) output feedback case, a possible choice 
would be to select corresponding elements from a full order solution 
as elements of an initial stabilizing reduced order gain matrix. Also, 
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the successive approximation technique used in DISLYP is similar to 
the one suggested by Kleinman (1974). The algorithm of this section 
is comparable to Hewer’s (1971) algorithm for solving the optimal dis- 
crete state regulator and is a discrete analog to the continuous al- 
gorithm developed by Levine and Athans (1970). 

4.2 Adaptive Sampling 

This section discusses adaptive sampling and gives a derivation 
for certain adapatlve sampling laws. As previously stated the pur- 
pose of adaptive sampling is to increase the sampling efficiency of a 
digital control system by varying the sampling rate with respect to a 
system parameter. In this regard the following definition is helpful. 

Definition . If is the value of a coat function associated 

with system A, and is the number of sampling instants used in 

generating that cost, then the sampling efficiency of system A, is 


Ha = 


1 

•Ja«a 


(4.2.1) 


Now the following design problem is stated. 

Problem Statement. Determine an algorithm that automatically ad- 
justs the sampling period (T^) based on some function of a scalar 
sampled signal y(t) or other system variables. That is, let 

Ti = fi(e(t)) (4.2.2) 


where 


e(t) = y(t) - y(ti) 


(4.2.3) 


and 


Ti = ti+i - ti 


(4.2.4) 


with the following constraint. 
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T < Tj ^ T 
'^tnin — — -^max 


(4.2.5) 


This constraint limits the variable sampling period, to an 

allowable range defined by the maximum processing rate of the digital 
controller (Tmin) and system stability requirements (Tmax) ♦ 

Problem Solution . The solution of this problem essentially fol- 
lows the approach of Hsia (1974) . Consider the cost functional 


Jl + Jj 


(4.2.6) 


where 


^ I Ue(t)|l^dt 


•^1 = 


(T.)‘ 


(4.2.7) 


and 

J2 = f2<Ti) (4.2.8) 

The cost J]^ is interpreted as the cost of incurring errors due to 
sampling, while J 2 is interpreted as the cost of taking samples. 

If the function sample per unit time, then 

for the region of constraint defined by equation (4.2.5), 

should be non-negative and monotonically decreasing. That is for 


'^ min — "ii — ”^max 


(4.2.5) 


then 

f2(Ti) >. 0 (4.2.9) 

and 

f2(Ti) > £ 2 ( 12 ). for Ti < T 2 (4.2.10) 

Also, to ensure a meaningful cost function a and b are restricted to 
the following sets of values. 
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a = {-1, 0, 1} (4.2.11) 

b = {1, 2} (4.2.12) 

Now an adaptive sampling law can be obtained by finding the that 

minimizes equation (4.2.6). 

To perform this minimization expand y(t) in a Taylor series about 

*=1- 

y(t) = y(ti) + y(t^)(t - t^) + . . . (4.2.13) 

Since 


’^max — ^ '’1 

and is chosen for stability of the system producing y(t) , it is 

assumed that higher order terms in the Taylor's expansion can be ne- 
glected. Then by (4.2.3) 


e(t) “ y(ti)(t - ti) 

Substituting (4.2.14) into (4.2.7) and integrating gives 

|y(ti)|^ 


j “ 


b + 1 


(Ti)^'^^''^ + f2(Ti) 


(4.2.14) 


(4.2.15) 


Differentiating this expression with respect to T^ yields 


3J 

3Ti 




(4.2.16) 


To determine a necessary condition for the minimization of (4.2.6) let 


•^=0 (4.2.17) 

3Ti 

From this necessary condition different adaptive sampling laws can be 
derived by choosing different functions for f2(T^). 
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For the adaptive sampling law used in this dissertation the fol- 
lowing choices were made 


a = 1 
b = 2 


and 


W - 


A(1 - BT^ + 




(4.2.19) 


Now made non-negative and monotonically decreasing by 


choosing A and B such that 

2 2 ^^ 

B Ti 

All - BT^ + — ^ I >. 0 


(4.2.20) 


and 


AB(BTj - 1) ^ 0 


for 


^min — *^1 — ^max 


The sampling period, T^, is constrained to be positive, 
(4,2.20) A is positive. From (4.2.21) B is positive and 

B < — 

-Ti 

Since the maximum value for Ti is T ^^ y > 


(4.2.21) 


(4.2,5) 
Then from 

(4.2.22) 


B = 


(4.2.23) 


max 


satisfies (4.2.22) and consequently £ 2 (^ 1 ) satisfies the conditions 
of (4.2.5), (4.2.9), and (4.2.10). 

Now from (4.2.16), (4.2.17), (4.2.18), and (4.2.23) 



80 


0 - |y(ti)p-|ii +5^4^^- - 

-* ^maxVmax 


(4.2.24) 


a[y(ti)) + 1 


(4.2.25) 


a 0 

AB^ 


(4.2.26) 


The adaptive sampling law is 


«[y(ti)] + 1 


^1 ^ Tjjln 


Tmtn» *^1 — "^min 


(4.2.27) 


The choice of a reflects a relative weighting between the cost of 
sampling and the cost of errors incurred due to sampling. As a de- 
creases the relative cost of sampling increases. Thus the number of 
sampling instants increases. 

As an example of a calculation for a let y(t) be approximated 
and normalized as 


. ^ ^ y(ti) - y(ti_i) 


(4.2.28) 


Assume that if yCt^.l) and y(t^) differ by at least 10% the sampling 
law should predict a minimum sampling period. From equation (4.2.25) 


max 

"i” - „ ^ + 1 

T? 

1 


(4.2.29) 
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Then 

a > (Tj. - 1)T? 100 (4.2.30) 

for all possible T^, where 


T = >> 1 


(4.2.31) 


Since the inequality of equation (4.2.30) must hold for all then 


a R; 100 


max 

^min 


It is remarked at this time that Smith (1971) has shown that the 


improvement in sampling efficiency may be highly dependent on the con~ 
tinuous process generating y(t), the presence of noise in the sampled 
signal, and the criterion defining performance. However, the adaptive 
control scheme described in the next section combines adaptive sampling 
with a control that is parametrically dependent on the sampling period. 
This minor sophistication over the systems studied by Smith gives im- 
provements in sampling efficiency that are not highly system or noise 


dependent . 

The selection of the system parameter used in the adaptive 
sampling scheme is highly system dependent. One possible choice is 
to use an output variable as the system sampling parameter that cor- 
responds to the slowest mode of the controlled system. Another pos- 
sible choice would be the average sampling error for all available 
system outputs. The large number of possible choices for system 
sampling parameters combined with the different adaptive sampling 
laws allow a great deal of design flexibility. 
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4.3 Adaptive Digital Control 

In the previous two sections both adaptive sampling and the dis- 
crete output regulator have been discussed. This section will combine 
these two control techniques into a simple adaptive digital control 
scheme. The objectives of design simplicity and practicality strongly 
influence the simplifications used to generate this control scheme. 

Essentially the control scheme Incorporates adaptive sampling 
and a proportional feedback gain matrix that is parametrically de- 
pendent on the sampling period. This is shown in block diagram form 
in Figure 4.1. The problem is to determine the functional relation- 
ship between the sample period, T, and the feedback gain matrix. 

Assume for the moment that the sample period is a constant. It 
is shown in Appendix D that if the continuous process is time-invariant 
the equivalent sampled-data system will be time-invariant. Assume also 
for the moment that full state feedback is available. Now the constant 
feedback matrix could be found by application of the steady-state dis- 
crete state regulator equations to the system in question. 

Now let the sampling rate vary with time. The sampled-data system 
is parametrically dependent on the sampling period and will be time- 
varying even though the continuous process is t ime- invar lant . If the 
sequence of sampling periods is known over the interval of control, 
the time evolution of the discrete system will be known. Therefore 
the feedback gain matrix could be generated by the time-varying dis- 
crete state regulator over the Interval of control. 

Now (still assuming full state feedback) Introduce adaptive 
sampling into the control picture. The sequence of sampling periods 
is now a function of a continuous system parameter. As such, the 
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Zero 

order 

hold 


Figure 4.1. - Adaptive control scheme. 


ORIGINAL PAGE IS 
OP POOR QUALOT 




84 


sequence will be unpredictable due to unknown system disturbances} 
changing environment, and noise < The most probable values of the 
sampling period will be the limit values T ma v« However, 

even the possible duration of a sequence of ^rain samples for in- 
stance, could not be determined. Thus the feedback matrix cannot 
realistically be determined by the time-varying regulator equations. 

The time-varying regulator solutions is parametrically dependent on 
the sample period and is computed backward in time. Thus sample 
periods that have not yet been predicted by the adaptive sampling law 
would be required for the solution of the time-varying regulator. 

To overcome this problem the following simplification will be 
made. For each sampling period predicted by the adaptive sampling law, 
the feedback matrix will be determined by the steady-state solution of 
the discrete regulator. Note that with this simplification either the 
full state regulator or the output regulator of the previous section 
can be used to determine the feedback gain, depending on the available 
outputs. Also, for a given linear, time-invariant continuous process, 
the simplification implies that there is a unique feedback matrix for 
each sampling period. Thus the feedback gains for the given process 
can be computed off-line for a sufficient number of sampling periods 
to establish a functional relationship between gains and sampling 
period. This relationship could easily be stored by a digital con- 
troller and used on-line to generate the control input to the con- 
tinuous process. This procedure would eliminate the prohibitive com- 
putation time needed to solve the discrete regulator problem on-line. 
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4»4 Example Problem 

This adaptive control scheme was applied to a linearized model of 
the Pratt and Whitney Aircraft F401 advanced technology two-spool turbo- 
fan engine. This model was Identified by Michael and Farrar (1973) . 

The F401 engine is being studied for use in a V/STOL (Vertical and short 
takeoff and landings) application and has variable exhaust, fan, and 
compressor geometries. The model selected for study here has a normal- 
ized fifth-order state vector and a normalized scalar control. The 
model represents the dynamical engine operation for small variations 
about a 73* power-lever angle setting. 

The fifth-order state represents turbine inlet temperature, com- 
bustor pressure, fan angular velocity, high pressure compressor angular 
velocity, and afterburner pressure, respectively. In this study the 
control represents jet exhaust nozzle area. 

In state-space form the model is 

x(t) * Fx(t) + Gu(t) 

where 


matrix is 
-34.0130 

-9.3030 

12.0370 

-2.3980 

-1.2540 

4.3890 

-38.7620 

-4.2210 

28.4800 

14.7290 

-4.7550 

2.2870 

-0.4000 

-1.5460 

-2.2000 

2.0460 

1.0620 

-0.7290 

-2.1500 

-0.6240 

4.1510 

-8.8140 

-0.1670 

7.4770 

1.0990 


and 
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The G matrix is 
0.7660 
0.0560 
0.1560 
-0.1370 
-4.7290 


The output vector is given as 
y(t) = Cx(t) 

where y(t) e R®. In this example each output element is one of the 
states so that each row of C has only one nonzero element and it is 
unity. 

The control objectives are assumed to be adequately described by 
the performance index 

J = -I [x’(t'lQx(t) + u'(t)Ru(t)]dt 

2 'Jo 


where Q » I and R = 1. To study the effect of adaptive sampling on 
the system's performance the adaptive sampling law 



Ti ^ Tjain 


T * < T • 


was implemented. The system parameter used to vary the sampling period 
was a combination of the engine outputs such that 


m 


e(t) = !yk(t) - yk<ti-i) 

k=l 
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where is the element of y(t) sampled at time tj^. Now 

the error derivative can be approximated as 

e(tj_) * e(ti)/Ti_i 

The limits on the sampling period were selected as 
^min “ 0.01 

and 


^max ^ ^ 


The Tjnij^ limit was selected as three times smaller than the real part 
of the smallest eigenvalue. The Tj„ax was selected to insure 

system stability. 

The engine model and the adaptive digital control scheme were simu- 
lated on an IBM 7094 digital computer. System performance was studied 
for the reduction of the initial condition disturbance 



to zero with different output configurations. In each experiment the 
performance index, J, the number of sampling Instants, N£, and the 
sampling efficiency 


^ “ Nf J 

were calculated. The results are summarized in Figure 4,2. 
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From these results it can be seen that sampling efficiency can be 
improved by a factor of 20 by the introduction of adaptive sampling with 
only a small increase in the system performance index. Note also that 
afterburner pressure is the most important of the states for control 
purposes. 


State 
available 
as output 

1,2, 3, 4,5 
1 , 2 , 3, 4 
1,2, 4,5 

2, 3, 4,5 

1 

2 

3 

4 

5 

5, 5 


Adaptive 

sampling 

parameter, 

a 


.0667 

.0667 

.0667 

.0667 

.0667 

.0667 

.0667 

.0667 

.0667 


Number of 

sampling 

instants, 

Nf 

500 

19 

24 

23 

14 
17 
12 
13 
13 

15 


Performance 

index, 

J 


.6711 

1.0415 

.6727 

.6733 

.9200 

.9148 

.9304 

.8847 

,6913 

,6760 


Sampling 

efficiency, 

1 


,00298 
. 05053 
.06194 
.06457 
.07764 
.06430 
,08957 
.08695 
.11127 
.09862 


State 




1 

2 

3 

4 

5 


Turbine inlet temperature 
Combustor pressure 
Fan angular velocity 

High pressure compressor angular , velocity 
Afterburner pressure 


Figure 4.2. - Adaptive digital control results for F401 
engine model. 



CHAPTER V 


APPLICATION OF ADAPTIVE DIGITAL CONTROL TO A JET ENGINE 

The purpose of this section is to describe the application of the 
previous results on identification and adaptive control to a jet engine. 
The second-order state-space operating point models given in Chap- 
ter III are used by the adaptive control scheme of Chapter IV to de- 
termine a piecewise constant control input for the digital simulation 
of a single spool turbojet engine. Thus, the adaptive nature of this 
engine control is twofold. First, the output feedback will vary when 
the engine model information reflects a change in the engine operating 
point. Second, the output feedback will change when the adaptive 
sampling law varies the sampling rate. 

The overall objective of this control is to provide rapid engine 
response to changes in the demanded steady-state operating condition 
while maintaining certain engine constraints. This engine regulation 
must be accomplished using only available engine outputs and with a 
minimum of computer control complexity and processing time. It has 
already been shown that the adaptive control scheme described in 
Chapter IV can satisfy the practicality constraints of the overall 
control objective. It remains to be shown that rapid response and 
engine constraint criteria can be satisfied by appropriate selection 
of the performance index and the output feedback variables. 
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5.1 The Engine Model 

From Chapter III the composite engine dynamics are written as 
6x =■ F(S)6x + G(S)6wf 

(5.1.1) 

6y = H(S)6x + E(S)6 w£ 

with 


y = y^ + i5y 


Wf “ Wfc 


(5.1.2) 


where S is the engine rotor speed and y^ and w^^ are the com- 
manded output vector and fuel flow, respectively. The commanded output, 
y^, is defined as the steady-state engine output corresponding to the 
constant engine fuel flow input, 

Recall that the state— space representation selected is second-order 
in the state and scalar in the control. The engine simulation can sup- 
ply as outputs temperatures, pressures, airflow, and thrust when given 
fuel flow, Wj> as input and the Mach number and altitude conditions. 

All results in this dissertation are given for sea-level, static con- 
ditions, i.e., the Mach number is zero and the altitude corresponds to 
sea-level. This condition is common among engine test studies. 

For an actual engine only certain physical outputs can be measured. 
Thus in this chapter only rotor speed, S, compressor discharge pres- 
sure, P^, and turbine inlet temperature, T-j., are assumed available from 
the engine simulation as possible feedback variables. Since the state 
is second-order, the output vector, y» is limited to first and second- 
order combinations of S, Pp, and Tj. Therefore, there exists six 
possible feedback combinations. Note that the elements of the com- 
manded output vector, y^., always correspond to the elements of the 


chosen output vector, y. 
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5.2 The Engine Control 

A block diagram of the control system and the engine is given in 
Figure 5.1. All the elements of this diagram have been discussed except 
the control gains, K^. Thus, the determination of these gains for cer- 
tain cost weightings is discussed. 

The Output Feedback Gains . The control objectives of the engine 
are assumed to be represented in the cost function 


J 


T r ^&y'Ct)Q6y(t) + R6w|(t)]dt 


(5.2.1) 


with R > 0 and Q = Q' > 0. This cost function and the state-space 
system of equation (5.1.1) can be equivalently rewritten as 


6x(t) = FjL6x(t) + Gj^dui(t) 

6yCt) = Hj^6x(t) + Ei6ui(t) 


(5.2.2) 


and 


J 


I j^*^^[6x’(t)Qi6x(t) + Ri6u|(t)]dt 


(5.2.3) 


where 

= F - G(E'QE + R)~^E’QH 

Hi = H - E(E'QE + R)~^E'QH (5.2.4) 

= H'(Q - QE(E'QE + R)“^ 'Q)H 
Rj^ = E'QE + R 

and 

6ui(t) = 6wf + (E'QE + R)”^E’QH6x(t) 


(5.2.5) 
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O - Variable rate sampler 
o - Summer 


n> 


- Schedule 


ISI 


- Miiltiplier 


Figure 5.1. - Engine ajid control system block diagram. 
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Since F, G, H, and E are functionally related to the engine rotor 
angular velocity, their counterparts in the equivalent system are also 
For the equivalent system of equations (5.2.2) to (5.2.5) and the 
sampling period T^, the discrete or sampled-data equivalent system is 
(see Appendix D) 

6x(k + 1) = A6x(k) + B6uj^(k) 

\D * ^ * 

6y(k) = C6x(k) + Dfiuj^Ck) 


Tii [6x '(k)Q'Sx(k) +26x’(k)M6ui(k) + 
^ k=0 


Defining an auxiliary variable 

6z(k) = C6x(k) (5,2 

and using the results obtained for the infinite time optimal discrete 
output regulator (t£ is assumed large with respect to system time 
constants) in Chapter IV, the output feedback control law becomes 


6ui(k) * -Kj^6z(k) (5.2.9) 

The piecewise constant control law over the continuous time interval 
tj^ < t < tj^+2 becomes 

6ui(k) = -Ki6z(ti) * -KiHi6x(ti) (5.2.10) 

Substituting the control transformation of equation (5.2.5) and the 
original system matrices into equation (5.2,10) and simplifying, the 
control law becomes 


Sw£(t) “ ~Kj5y(ti) 


( 5 . 2 . 11 ) 


where 
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* (Ki + R~^E'Q)(I - EKi)"^ (5.2.12) 

Therefore, the gain matrix, Kj^, is both a function of the sampling 
period, T^, and, implicitly, the engine rotor speed, S. It is remarked 
that the feedback matrix can be obtained by either scheduling system 
matrices as a function of S and computing the feedback matrix at each 
sampling instant on-line, or by scheduling the feedback matrix itself 
as a function of S and (The feedback matrix would be computed 

of f-line for sufficiently large intervals of S and and stored.) 

The first technique would require less storage than the second but 
would require more on-line computing time. The second technique 
would require less on-line computing time but more storage capacity 
than the first. Since the critical element in most digital computer 
control systems is on-line computing time, the second technique is 
preferable. Since a digital simulation rather than a real time en- 
gine is the process to be controlled in this paper, on-line computing 
time is not a problem. Thus, the first technique Is used in this dis- 
sertation for the resultant savings in computer storage requirements. 

The Cost Weighting Matrices . The cost function given in equa- 
tion (5.2.1) represents a quadratic weighting of the output and control 
energy excursions from the steady-state operating line. The choices of 
Q and R, the weighting matrices, will reflect the relative importance 
of Incurring output and control errors. Since each element of 6y(t) 
and 6wf(t) vary in relative magnitudes, scaling was introduced to 
facilitate the choice of the weighting matrices. 

To each sampled engine rotor velocity, N = there corresponds 
a steady-state operating condition of the engine. The value of the 



95 


output vector and control at this condition are yjj and Wfjj. If a 
matrix Qjj is defined as 

On = DIAGtyj,] 

the percentage change in the output error vector can be written as 


6ypc(t) = Qii^SyCt) ^5.2.1 

Similarly the percentage change in the control error can be written 

SWfpcCt) = 6Wf(t)/wjjj (5.1.1: 

In terms of the percentage changes in error the cost function of equa- 


tion (5.2.1) becomes 


1 

■ 2 Jo 


[6yj;^(t)Q6yp^(t) + Rfiw^pjdt 


(5.2.16) 


Now the choice of Q and R reflect relative weightings for the per- 
centage changes in output and control errors. The original weighting 
matrices, Q and R, are determined by 


Q = Qn^<?Qn^ 
R = Rj^^RRjj^ 


(5.2.17) 


5.3 The Adaptive Sampling Law 

From the previous chapter the adaptive sampling law was defined as 


r 


■■max 


T. > T, 


^1 „ Jae(ti) + 1 


min 


(5.3,1) 


Vjmln’ 


^i — ^min 

where e(t^) is the difference between the past and present sampled 


values of some system parameter. 
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For the application of this adaptive sampling law to the engine 
controller, engine rotor velocity was selected as the system parameter 
from which the sampling period is predicted. This choice was made for 
two reasons. First, engine rotor velocity is already required by the 
feedback gain schedule and is therefore already available. Second, 

S, as the slowest and smoothest responding engine variable, gives the 
best indication of the "dynamic position" of the engine relative to 
steady-state. 

It is remarked at this point that if a system is in the steady- 
state for a long period and a sudden disturbance occurs, the slowest 
system output might not decrease the sampling period fast enough to 
handle the disturbance. This problem can be handled by conservatively 
choosing the T ^^v limit or by using both the slow and fast system 
outputs in some weighted combination to predict T^. Such problems 
are not considered here. 

Approximating the scaled or percentage error derivative as 


®pc^^i^ 


S(tj^) - S(ti_i) 

Ti-lS(ti-i) 


( 5 . 3 . 2 ) 


the adaptive sampling law is 


r 


*^max 


^1 


a(4p^(ti»^ + 1 


Tl > Tj„in 




min* 


( 5 . 3 . 3 ) 


T-r < T, 


mxn 


Preliminary engine simulations established the sampling limits as 


Tmln = 0-001 


*^max 


0.025 


( 5 . 3 . 5 ) 



97 


and the weighting term as 
a = 4x10^ 


(5.3.6) 


5.4 The Simulation Results 

The first simulation results were obtained for the case of rotor 
velocity as the only variable for feedback. In this case y = S. The 
control task was to accelerate the engine from steady-state at 90% to 
steady-state at 104.5% design speed (100% is 36 960 rev/mln) . The 
commanded fuel flow which represents this change in the steady-state 
engine condition is 


w 


fc 


To. 09244 Ibm/sec 
1.0.15 Ibm/sec 


0 < t < 0.1' 


0.1 < t .< 2 


(5.4.1) 


As a means of comparison, baseline results with a constant sampling 
period of T = 0.001 seconds were simulated with ~ ^ sec. The 
scaled weighting matrices were selected as 

Q « 10 

(5.4.2) 

R = 1 


The respective trajectories for rotor speed, S, compressor discharge 
pressure, P^, turbine inlet temperature, Tj', engine thrust, F^, and 
fuel flow, Wf, are given in Figures 5.2 to 5.6 for the commanded fuel 
flow Wfj.. For future reference label this simulation test "Case 1." 

For the baseline (Case 1) and subsequent simulations the weighted 
cost function 


= Y ^ [sy’(t)Q6y(t) + R6w|(t)]dt 


(5.2.1) 


and the unweighted coat function 
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are calculated* Using the definition of sampling efficiency from Chap- 
ter IV, the weighted and unweighted sampling efficiencies 


riw 


Vf 




(5.4.4) 


(5.4.5) 


where Nf is the number of sampling instants are calculated for each 
simulation. These indices are summarized for each simulation at the 
end of this chapter in Figure 5.27. 

Next, adaptive sampling was added and the simulation repeated for 
the commanded fuel flow of equation (5.4.1). Call this "Case 2." 
Again, 


Q * 10 
R = 1 


(5.4.6) 


The resultant trajectories are plotted with the previous baseline re- 
sults to show the effect of including adaptive sampling in Figures 
5.7 to 5.11. The results show no visible differences except in the 
turbine inlet temperature plot. A comparison of the cost and efficiency 
indices shows approximately a 10% Increase in the cost and a 100% im- 
provement in sampling efficiency. 

To examine the effect of a change in the cost function, the cost 
function weighting was changed to 


Q = 50 
R = 1 


(5.4.7) 
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Figure 5*7^ - Case 2 and Case 1 engine accelerations. 
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Figure 5.11. - Case 2 and Case 1 engine accelerations. 
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and the simulation repeated with adaptive sampling and the commanded 
fuel flow of equation (5.4.1). Label this simulation "Case 3." These 
results are plotted with those of Case 2 and are given in Figures 
5.12 to 5.16. From the figures it is clear that the Case 3 control, 
as expected, accelerates the engine faster than the Case 2 control. 

For example, rotor speed is within 1% of its final steady-state value 
in 0.6 sec for case 3 while in Case 2 it required 1.1 sec. For en- 
gine thrust a similar analysis gives acceleration rise times of 
0.4 sec for Case 3 and 1.2 sec for Case 2. The penalty for this im- 
proved acceleration of Case 3 is higher turbine inlet temperatures and 
additional fuel flow requirements, A comparison of the cost and effi- 
ciency indices shows that the weighted and unweighted costs are higher 
and the sampling efficiencies lower for Case 3 than for Case 2. 

From the trajectories already presented, it is seen that the 
adaptive control described in this paper can rapidly accelerate the 
engine from one operating point to another. However, in realistic en- 
gines the temperature history of the turbine blades is also of utmost 
importance. High turbine inlet temperatures, Tt, can cause short life- 
times or even outright failures of the turbine blades. Thus another 
important control consideration is the limiting of turbine inlet tem- 
perature. To evaluate the ability of the adaptive digital control 
scheme to limit turbine temperature and still effectively accelerate 
the engine, a simulation test, called "Case 4," is devised as follows. 

Suppose that the temperature limit for the single spool turbojet 
engine is 1900® R. The control configuration must now accelerate the 
engine from 90% to 104.5% design speed without violating the 1900® R 
temperature constraint. The commanded fuel flow is again as in 
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Figure 5-1^. - Case 3 and Case 2 engine accelerations. 





Figure 5 . 15 * ~ Case 3 and case 2 engine accelerations. 
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Figure 5*l6. - Case 3 and Case 2 engine accelerations. 
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equation (5.4.1). This temperature constrained acceleration is accom- 
plished by initially controlling the engine as in Case 2 with engine 
rotor speed feedback, i.e., 

y = S (5.4.8) 

and 

Q «= 10 

(5.4.9) 

R “ 1 

until the 1900® R level is exceeded. Then the control configuration is 
changed to Include both rotor speed, S, and turbine temperature, Tj, 
feedback. The output vector is now 



and the cost weighting matrices are chosen as 



R * 1 


(5.4.10) 


(5.4.11) 


As before this weighting still penalizes the rotor velocity error, but 
now it also penalizes large temperature errors. In addition it re- 
flects a greater importance on the conservation of control energy. 

Once the temperature constraint is satisfied, the control configuration 
reverts to the original (Case 2) feedback arrangement. The results 
for Case 4 are plotted along with those of Case 2 in Figures 5.17 
to 5.21. From these figures it is seen that the turbine temperature 
is limited to 1900“ R with only a small penalty in acceleration when 
compared to the Case 2 acceleration. A comparison of cost and 



Time (sec) 


Figure 5«1T* - Case and Case 2 engine acceleration. 
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Figure 5.19. - Case k and Case 2 engine accelerations. 
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Figure 5-20. - Case 1| and Case 2 engine accelerations. 
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efficiency indices for Case 2 and Case 4 shows only minor variations in 
these values. 

A final simulation test. Case 5, is now developed to determine if 
the addition of a second feedback variable, P^., can improve engine ac- 
celeration, Again the commanded fuel flow input is given by equa- 
tion (5.4.1), The feedback variables are 

(5.4.12) 

and the weighting matrices are selected as 



R = 1 

The result of Case 5 are plotted along with the trajectories of Case 2 
in Figures 5.22 to 5.26. 

These figures show that no appreciable improvement in accelera- 
tion is obtained by the addition of the feedback variable Turbine 

inlet temperature is lowered somewhat, but fuel flow is increased in 
Case 5 when compared to Case 2. In addition, a comparison of the cost 
Indices shows a sizeable increase in the cost from Case 5 to Case 2 
while the sampling efficiency went down. However, it should be noted 
that the weighting selected does not emphasize the elimination of 
speed error as heavily as previous test simulations. In this regard 
it should be mentioned that the weighting matrices selected in each 
case were probably not the best possible. If required, more time 
could be devoted to achieving better results by additional test simu- 



lations . 





-Turbine inlet temperature (^R) 
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Figure 5*2^. - Case 5 sind Case 2 engine accelerations. 


-Engine thrust (lb) 


118 



Time (sec) 


Figure 5 * 25 , - Case 5 and Case 2 engine accelerations. 
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Figure 5«26. - Case 5 and Case 2 engine accelerations. 



ORIGINAL PAGE iS 
OF POOR QUALITY 


Simulation 

test 

Case 1 
Case 2 
Case 3 
Case 4 


Case 5 


Output 

feedback 

variables, 

y 

Weighting 

matrices, 

Q,R 

Number of 

sampling 

instants, 

% 

, s 

^10,fcl 

1928 

s 

Q=10,^l 

878 

s 

^50,^1 

621 

s 




Unweighted 

cost, 

^u 

Weighted 

cost, 

"^w 

Unweighted 

sampling 

efficiency, 

Weighted 

sampling 

efficiency 



.246E-1 

. 315E-2 

.211E-1 

.165 

.278E-1 

.337E-2 

. 409E-1 

.338 

.591E-1 

.256E-1 

.273E-1 

. 628E-1 


or 




.297E-1 .279E-2 .374E-1 .398 


.313E-1 .658E-2 .380E-1 .180 


Figure 5.27. - Summary of simulation results. 
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5.5 The General Control Procedure 

If this control procedure were to be applied to a different engine > 
the following steps would be required. First, identify operating point 
models from data generated by the new engine at enough points to ade- 
quately describe the engine dynamics. Second, select the appropriate 
outputs and schedule the commanded output against the commanded input. 
Third, apply the adaptive control scheme of Chapter IV to the iden- 
tified model dynamics for appropriately selected cost functions, 
sampling periods, and operating points. Fourth, schedule the result- 
ant feedback gains as a function of operating point and sampling 
period and store this function in the control computer. Finally, se- 
lect a system parameter to be used by the adaptive sampling law to 
predict the sampling period. The total control system is now as shown 
in Figure 5.1. 



CHAPTER VI 


SUMMARY OP RESULTS AND CONCLUSIONS 
In Chapter I the need for new concepts in the control of jet en- 
gines was discussed. It was pointed out that modern control theory 
had significant advantages in the development of these concepts. 
Specifically these advantages were the ability of modern control 
theory to design multivariable control systems that take advantage of 
loop interactions and the systematic way in which the control systems 
were designed. Specific research objectives were the identification 
of dynamic engine models from realistic engine data with a minimum 
of a priori assumptions, and the development of computer control al- 
gorithms that were both efficient and practical. 

In Chapter II a brief description of the physical characteristics 
of air-breathing gas turbine engines was given. A summary of the 
basic concepts of engine control was also presented. 

In Chapter III the identification of a low-order dynamic state- 
space model for the single spool turbojet engine was described. It 
was shown that a technique developed by Tse and Welnert could identify 
steady-state linearized operating point models with a minimum of 
a priori assumptions. For this technique both the model order and 
model parameters were determined in a noniterative fashion from re- 
alistic data generated by a digital computer dynamic simulation. 
Gradient techniques were used to complete the identification of the 
model and a comparison of the composite model and the engine 
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simulation was performed. 

In Chapter IV discrete time algorithms were developed for effi- 
cient and practical computer control of linear systems. The optimal 
discrete output regulator problem on a semi— inf inite interval was pro- 
posed and the necessary conditions for optimality derived by Lagrangian 
techniques. An algorithm for the solution of these necessary condi- 
tions was presented along with the computer listings of the required 
programs. Also, adaptive sampling and its ability to improve sampling 
efficiency were discussed. Next, adaptive sampling and the optimal 
discrete output regulator were combined to form an adaptive digital 
control scheme that was applied to a fifth-order linearized engine 
model. Both the sampling efficiency and control degradation under 
different output feedback configurations were studied. 

In Chapter V the adaptive digital control scheme was applied to 
the operating point models of Chapter III and the resultant adaptive 
configuration used to control the single spool turbojet engine simula- 
tion. The twofold adaptive nature of the feedback matrix as a func- 
tion of rotor angular velocity and sampling period was described. 

The results of the simulations showed that rapid engine acceleration 
from one operating point to another could be achieved with this 
adaptive control scheme using only rotor velocity feedback. The ef- 
fect of a change in the weighting matrices on the acceleration time 
was studied along with the improvement in sampling efficiency. It 
was also shown that the addition of temperature feedback could be 
used to limit the maximum turbine inlet temperature. Finally, the 
control configuration of rotor velocity and compressor discharge pres- 
sure was simulated to study the effect of an additional feedback 
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variable on acceleration time and the regulation of the pressure vari- 
able. The weighted and unweighted costs and sampling efficiencies for 
each simulation were summarized in Figure 5.27. 

From this study it can be concluded that modern control theory 
can be successfully applied to jet engine control. In particular using 
linearized operating point models to describe engine dynamics, the 
adaptive digital control scheme of Chapter IV can successfully control 
a jet engine using available outputs in a computationally efficient 
manner . 

Specific achievements of this dissertation include 

1. The identification and verification of a second-order state- 
space model for a turbojet engine using an identification method by 
Tse and Weiner t and realistic engine data. 

2. The derivation of the necessary conditions of optimality for 
the optimal discrete output regulator with crossweighting in the per- 
formance index on the semi-inf inite time interval. 

3. The development of a computer algorithm to solve the neces- 
sary conditions in (2). 

4. The combination of adaptive sampling and the optimal discrete 
output regulator into an adaptive control scheme. 

5. The application of the adaptive control scheme of (4) and the 
model of (1) to the control of a jet engine. 

There are several worthwhile extensions to this research. One 
is the development of an on-line technique that identifies engine dy- 
namics. Such a technique would minimize the time required for the 
initial control system design of a number of engines with the same 
configuration. In this regard the Tse and Weinert Identification 
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technique may be used to initialize an on-line technique. A second 
extension is to increase the number of control variables (to include 
exhaust nozzle area, e.g.) and the flight conditions for which the 
control is designed. Since an engine operates over a range of flight 
conditions, a realistic control must be designed for all these condi- 
tions. Third, different adaptive sampling laws could be applied to 
the adaptive digital engine control scheme to determine if better 
sampling efficiencies can be obtained. Finally, the benefits of using 
additional sensed outputs as feedback variables could be compared to 
the cost of additional sensors by the output regulator formulation of 
Chapter IV. Similarly, the effect of sensor failures on engine per- 
formance could be evaluated by the output regulator formulation. 



APPENDIX A 


TIME SERIES ANALYSIS 

Time series analysis is used to determine the parameter set 
{pi, described in the section on the Tse and Weinert identifica- 

tion method. The following derivation from Tse and Weinert (1973) 
shows how this is accomplished. Let 

R(a) = E{y(k + o)y'(k)}, o = 0, 1, 2, . . . (A.l) 

E = E{z(k)z'(k)} (A. 2) 

Given the system 

z(k + 1) “ Az(k) + Bv(k) 

(A. 3) 

y(k) = Cz(k) + v(k) 

where z(k) eR", y(k) eR™, and v(k) is a zero mean Gaussian noise 
process with covariance 

E{v(k)v'(j)} “ Q6kj (A. 4) 

The following equations can be derived 

E = AEA* + BQB’ (A. 5) 

R(o) = CE{z(k + a)z'(k)}C' + CE{z(k + 0)v’(k)}, a > 0 (A.6) 

But 

E{z(k + o)z'(k)} = A^Z (A. 7) 

and 

E{z(k + o)v'(k)} “ A'^'^BQ, o > 0 (A. 8) 


126 



127 


R(a) = CA'^-1(AEC + BQ) - CA®"^S 


(A. 9) 


S = Arc' + BQ 


R(0) « CEC’ + Q 


Let rij(o) be the element of R(a) , and Sj be the j 


column of S. Then 

. V t . o-l„ 

rij(a) - C.A s^, 


a > 0 


where is the row of C. 

Assuming that A and C are in the canonical forms given in the 
chapter on identification, the parameters {pi> Bijk) are related by 


c/i - C 

j=l k=0 


eijkcjA^ 


if Pi > 0 


1-1 pj-i 

ci = iz: eijkcjAk 

j=l k=0 


if Pi = 0 


Now using equations (A. 13) and (A. 14) 

^ P&-1 

rij(Pi + t) = / , ) i Bijjij t^A^A’^'^Sj , 

^ A =1 k =0 


P^ > 0 


i-1 1^-1 

rij(pi + t) “ ZZIZZ! ’ 

£=1 k »0 


P-i = 0 


where t = 1, 2, • ^ and 


rijCPi + t) = CiA^^A'^"^Sj 



Then using (A. 12) 


i Pj^-1 

ry (Pi + t) -Cn; (k + t) . Pi > 0 (A.18) 

S,=l k=0 ^ 

i-1 P£-l 

*■ 1-1 ^Pi + t) = 3 i + t), Pi = 0 (A. 19) 

£«1 k=0 


Recall that the parameter set is determined from the identif labil- 

ity matrix. The identif lability matrix is composed of elements of R(o). 
Once the parameter set is determined, equations (A.18) and (A. 19) can 
be solved for{^j]^}. Thus the parameter set {P^, is determined 

from the matrix R(a) . This matrix can be estimated by time series 


analysis as 



APPENDIX B 


LOGARITHM OF A MATRIX 
Consider the matrix equation 

eX = A (B.l) 

All the solutions to this equation are called logarithms of A and 
are denoted by In A. 

The characteristic values Xj of A and Cj of X satisfy 

Xj = (B,2) 

Assume that det{A) 0 and X^ X j , that is the eigenvalues are 

distinct. Let 

A = PAP"^ (B.3) 

where A is a diagonal matrix and P is a similarity transformation 
matrix. 

A = DIAG[ajj] = DIAGlXj] (B.4) 

Then from Gantmacher (1959) the matrix X satisfies the equation 

X = P In AP~^ (B.5) 

and the logarithm of A is 

In A = DIAG[ln(Xj)] « DIAG[Cj] vB.6) 
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APPENDIX C 


COMPUTER SIMULATION SUBPROGRAMS FOR THE SOLUTION OF THE 
DISCRETE OPTIMAL OUTPUT REGULATOR 
CLSDLP 
DISLYP 
DITORF 
MULT 
RICATT 
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C SUBROUTINE MULT 

C 

C PURPOSE 

C TO COMPUTE PRODUCT OF TWO MATRICES 

C 

C DESCRIPTION OF PARAMETERS 

C ALPHA- N X L REAL MATRIX 

C BETA - L X M REAL MATRIX 

C GAMMA- N X M REAL MATRIX 

C N - NUMBER OF ROMS IN ALPHA 

c M - number of columns in beta 

c L - number of columns in alphacrows in beta> 

c 

subroutine multi alpha, bet a. gamma, N fM,L» 

DIMENSION ALPHA(11,11),BETA(11,U) fGAMMAI H ,11) 

DO 10 I*1*N 

DO 10 J<1,M 

GAMMA(I,J)=0»0 
DO 10 K=1.L 

10 GAMMA! I,J)-GAMNA( I, JM-ALPHA! I ,K) *BETA (K , J ) 

RETURN 

END 


C SUBROUTINE CLSDLP 
C 

C PURPOSE 

C TO CALCULATE THE CLOSED LOOP MATRIX 

C AA AND ITS TRANSPOSE 

C AA * A - B*F»C 

C 

SUBROUTINE C LSDLP < A ,B.F, C. AA . AAT ,N*M,L,NMAX > 

DIMENSION A!NMAX,1) ,B(NMAX,1) ,F(NMAX,1) ,C (NMAX, I > , AA( NMAX , I ) 
DIMENSION AAT(NMAX,1) 

CALL MULTIF.C.AAT.M.N.L) 

CALL MULT(B, AAT,AA,N,N,H) 

DO 10 1=1, N 

DO 10 J=1,N 

AA< I,J>= A(I ,J> - AA(I,J) 

10 AATCJ,n = AA< I,J) 

RETURN 

END 
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SUBROUTINE DISLTP 
PURPOSE 

TO SOLVE THE DISCRETE LYAPUNOV EQUATION 
P s A*P*A» + S 

METHOD 

SUCCESSIVE SUBSTITUTION 
SUBROUTINE D I SL YP ( P . A , S , ERROR t Nf NMAX) 

DIMENSION P(NMAX,l>»A(NMAX.l).S(NMAX.l) .F ( 1 1 . 11 ) . ATI U. 11 ) 
DIMENSION Aim.ll) .smi.ii) 

DO 1 1=1, N 

00 1 J=1.N 
A1(I,J)=A(I, J) 

1 S1(I,J)*S( I, J) 

1=0 

5 SUM=0, 

1 = 1*1 

00 10 1=1. N 
DO 10 J=1,N 
P( I • J> = Sl( I,J> 

10 ATIJ.I) = Al(I.J) 

CALL MULTU1,P,$1,N,N,N) 

CALL MULTI Sl.AT.F.N.N.Nl 
00 20 1=1, N 
DO 20 J=1,N 

20 Sltl.J) = PI l.J) ♦ FII.J) 

C CHECK FOR CONVERGENCE 
DO 30 1 = 1, N 
00 30 J=1,N 
F< I,J) = All I,J> 

30 SUM = ABS(Sin,J) - P(I,J)> * SUM 
IF(SUM.LT*ERROR> RETURN 
CALL MULT(F,F,A1,N,N,N) 

IF< I .LT.25> GO TO 5 

RETURN 

END 


OF POOK ClUALlTY 
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SUBROUTINE OITORF 
PURPOSE 

TO SOLVE the following EQUATION 

F * INVIR ♦ B**P*B)*IB**P’t'A*H')*EL*C**INVIC*EL*C' > 

FOR THE DISCRETE. INFINITE TIME, OUTPUT FEEDBACK 
REGULATOR MATRIX F 

SUBROUTINE 0 ITORF ( R ,8 ,P, A,EL, C i AM,F ,N ,M,L ,NMAX > 

DIMENSION R(NMAXfl) ^BINMAXyll tPINMAXtU ,A(NMAX« 1) 

DIMENSION EL(NMAX*U.C(NNAX.1).AM(NMAX« 1> 

DIMENSION BTKIll,ll),ELCT(ll,ll),DUMHllrU)»BTUl,ll) ,CTUl,ll> 
DIMENSION LWORK(ll) iMWORKUl) 

KMAX*11 
DO 10 1=1, N 
00 20 J=1»M 
20 BT<J. I)»BU. J) 

DO 30 K«1,L 
30 ctu,k»=c<k. n 
10 CONTINUE 

CALL MULTIEL,CT,ELCT,N,L,N) 

CALL MULT<8T.P.BTK.M«N«N> 

CALL MULT(C,£LCT,0UM1,L,L,N) 

C INVERT THE MATRIX OUMls C » EL ♦ C* 

CALL ARRAY(2.L.L.KMAX.KMAX.DUM1,DUM1) 

CALL MINV(0UM1,L,0ETC,LW0RK,MWQRK) 

CALL ARRAY (lfL,L,KMAX,KMAX,DUMl,DUMl» 

CALL MULTCELCT.OUMl.CT.N.L.LI 
CALL MULT(BTK,A,DUM1,M,N,NI 
DO 50 1=1, M 
DO 50 J=1,N 

50 0UM1(I,J) = OUMKItJ) ♦ AMU, I) 

CALL MULTIDUMl.CT.ELCT.M.L.N) 

CALL MULT<8TK,B,0UM1,M,M,N) 

DO 40 1=1, M 
DO 40 J=1.M 

40 BT( I,J)=R( I, J>+DUM1(I,J) 

C INVERT THE MATRIX BT » R ♦ 8**K*B 

CALL ARRAYCZ.M.M.KMAX.KMAX.BT.8TJ 
CALL MINVI 8T ,M, DETBT,LWORK,MWORK> 

CALL ARRAY(l,M,M,KMAX,KMAX,BT,8TJ 
CALL MULTtBT.ELCT.F.M.L.Ml 
RETURN 
END 
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C SUBROUTINE R ICATT 

C 

C PURPOSE 

C TO SOLVE THE OUTPUT FEEDBACK REGULATOR 

C PROBLEM FOR THE STEADY STATE DISCRETE PROBLEM USING 

C the discrete analogue of ATHANS and LEVINE'S METHOD 

C 

subroutine RICATTtA.B.C >Q*AM. R.F.N«M*L»NMAX> 

DIMENSION A(NMAX,1>,8INHAX*1) fCINMAX, 1),Q(NMAX,X>,RINMAX,1> 
DIMENSION AM(NMAX,1» »FINMAX,l> 

DIMENSION P( 11«11)*DDC11 *11>*AAU1>LU«AAT(11.11>.QI0ENT(11»11) 
DIMENSION el (ll,U>,OEm,ll)»DFlllfll> ,DG«ll,ll) 

DATA NMAAtNMF,NML»NMP /2HAA, IHF, IHL.IHP / 

ERROR * FL0ATCN)»»2»l»E-5 
ERR = FLOAT! MI^FLOATILI^I^E-S 
DO 5 1=1. N 

5 QIOENT( I,I »»1.0 
100 CONTINUE 
C 

C CALCULATE THE CLOSED LOOP SYSTEM MATRIX AA= A - B*F*C 
C 

CALL CLSOLPt A.B.F.C.AA.AAT.N.M.L.NMAXI 

C 

C GIVEN A VALUE FOR F CALCULATE THE P MATRIX 

C 

CALL MULT(F,C»0D,M,N»L) 

DO 30 1=1, M 
DO 30 J=1.N 
30 P( J, I J=DD(I, JI 

CALL MULTIP.OO.OF.M.N.M) 

CALL MULTIP, DF.DG,N,N,M) 

CALL MULT! AM,DO»P,N,N,M> 

DO 40 1=1, N 
DO 40 J=1,N 

40 DO! I,J >»0G( 1 tJI-i-Od , J>-P! I, J)-P! J, I) 

CALL DISLYPi P.AAT.DD, ERROR, M.NMAX) 

C 

C COMPARE PAST AND PRESENT ALGORITHM SOLUTIONS 
C 

SUM = 0, 

DO 10 1=1, N 
DO 10 J=l.N 

HOLD = OE! I, J> - PI I,J» 

IFIP! I.J) .NE.O. » HOLD = HOLD/P<I,J» 

SUM = ABS(HOLD» ♦ SUM 
10 OE! I , J)=P< I, J) 

IF(SUM#LT,10,*EP.R0R» GO TO 200 
50 CONTINUE 
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C 

C CALCULATE THE L MATRIX 

C 

CALL CLSDLPt A*B.F.C.AA.AAT.N.M.L.NMAX) 

CALL OISLYP(EL,AAtQIOENT,ERROR,N,NMAX) 

C 

C CALCULATE A NEW F MATRIX CALLED DO 

C 

CALL DlTORFIR,B»PfA»EL*C» AM,00,N,M,L.NMAX> 

C 

C COMPARE PAST AND PRESENT GAIN MATRICES 

C 

SUM= 0*0 
DO 20 1*1, H 

DO 20 4=1. L 

IF(FtI.4>.NE.O* > HOLD = HOLD/F(I,J> 

HOLD = DD(I,J> - FCI,J) 

SUM » A6SCH0LD) ♦ SUM 
20 F( I,J}=OD( It J) 

IF(SUM*6T.ERP> GO TO 50 
GO TO 100 
200 CONTINUE 
RETURN 
END 
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APPENDIX D 


SAMPLED-DATA SYSTEMS 

Many processes are naturally modeled as continuous time processes. 
However, the introduction of a computer as the principal control ele- 
ment requires that time be quantized to match computer processing time. 
This combination of continuous process, digital controller, and sampling 
operation between process output and controller input is denoted as a 
sampled-data system (Levis, et al., 1971). A sampled-data system is 
given in Figure D.l. 

To represent a sampled-data system mathematically let the continu- 
ous process be described as a time- invariant linear vector differential 
system of equations 

X = Fx + Gu 

(D.l) 

y = Hx 

where x e R“, y ^R®. and u e R9. Since u(t) is the output of a 
digital controller, a sample and hold constraint is imposed. on u(t). 
Therefore, the control is assumed piecewise-constant with changes in 
its value only at sampling Instants tj^. Thus 

u(t) = u(t£> = uj, for ti :< t ± t^+1 (D.2) 

The continuous system of (D.l) is now discretized as follows. 

First the trajectory x(t) for the system of equations (D.l) is given 
by 
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Zero 

order 


hold 


Figure D-1. - Saunpled-data system. 





t 


(D.3) 


x(t) 


/(t-ti) 


x(ti) + 
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gF(t-T)gy^ dT 


Define 

T * - ti 

A =“ A(T) « 

pT 

B = B(T) = / e^^G dt 
'-'O 


(D.4) 

(D.5) 

(D-6) 


H “ C CD. 7) 

x(ti^.j^) * x(i + 1), x(tj^) = x(i), y(tj^) = y(i) (D.8) 

Then the discrete equivalent of (D.l) is 


x(i + 1) = A(i)x(i) + B(i)u(i) 

(D.9) 

y(l) = Cx(i) 

If the sample rate, T, Is constant for all i, then the system (D.9) 
becomes time- Invariant . 

Given a quadratic cost function in terms of the continuous vectors 
y(t) and u(t) 


J - [y'Ct)Qy(t) + u’(t)Ru(t)]dt + Y y’(tf)QjyCtf) (D.IO) 

Jto 


where 


>. 0, Q i 0, and R > 0 


(D.ll) 


Also, 


Qf “ Qf. Q’ “ Q. and R' = R 


(D.12) 


Finally 
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Qf + Q 0 


(D.13) 

The discrete equivalent is 



N-1 

J = x'(N)Qjjx(N) + -r ) (x*(k) ^x(k) 

^ k=0 

+ 2x*(k)Mu(k) 



+ u’(k)Ru(k)} 

(D.14) 

where i = 0, 1, . . N - 1, when 



T “ H+i - ti 


(D.4) 

y(k) » Cx(k) 


(D.9) 

Qn = CQfC' 


(D.15) 

/oT 

Q(T) = / A'(t)CQC’A(t)dt 

do 


(D.16) 

M(T) = r A' (t)CQC'B(t)dt 

Jo 


(D.17) 

and finally 



R(t) ■= TR + P B'(t)CQC'B(t)dt 

Jo 


(D.18) 


The weighting matrices Qj|, Q, R, and M are parametrically de- 
pendent on the sampling period. Thus if the sampling period varies, 
the discrete system would be time-varying even though the continuous 
system is time-invariant. 

The matrix A(T) , being a fundamental matrix, is nonsingular. 
Since Q(R) is symmetric and positive semi -definite (definite) , it is 
easily shown that Q(R) is also symmetric and positive semi-definite 
(definite) . 
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